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In this review, I present the description of the early stages of heavy ion collisions 
at high energy in the Color Glass Condensate framework, from the pre-collision 
high energy nuclear wavefunction to the point where hydrodynamics may start 
becoming applicable. 


1. Introduction 

Heavy ion collisions pose a challenge for Quantum Chromodynamics (QCD) because 
a comprehensive description of these collisions involves a mix of hard short distance 
phenomena, and non perturbative long distance soft physics. Even the aspects of 
these collisions where a hard scale may justify a weak coupling approach are not 
perturbative in the naive sense of a strict loop expansion. Indeed, resummations 
are often required even though the coupling is small, usually because the bosonic 
constituents of the system have a large occupation number that may compensate 
the smallness of the coupling. 

One of the areas where a weak coupling QCD approach is expected to be most 
effective is the description of the early stages of a heavy ion collision. The term 
“early stages” usually encompasses the description of the relevant degrees of freedom 
in the wavefunctions of the two projectiles prior to their collision, the interactions 
that happen during the very brief duration of the collision itself, and the subsequent 
evolution of the produced gluons and quarks shortly after the collision. Roughly 
speaking, the temperature (or the fourth root of the energy density if the system 
is not yet thermalized) of the system can serve as a measure of the applicability 
of weak coupling techniques, since this scale sets the value of the running coupling 
constant. 

A lot of progress has been made in the last 20 years in understanding how to 
apply QCD to these collisions. The starting point was the realization that a hard 
scale naturally emerges from the non-linear interactions among the gluons when 
their density is large, 1-3 as is the case in the wavefunction of high energy hadrons 
or nuclei. Thus, the bulk of particle production in high energy heavy ion collisions 
is amenable to a weak coupling treatment. The formalism for doing this -known 
as the Color Glass Condensate (CGC)- was progressively established and refined 
during this period, and has by now reached a mature state allowing quantitative 
and systematic calculations. 
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An outstanding problem, that has not yet reached a satisfactory state of under¬ 
standing, is the transition from the Color Glass Condensate description to a more 
macroscopic description such as hydrodynamics. The main question is to explain, 
within the CGC framework, why an hydrodynamical description is possible in the 
first place. In other words, why does the system produced in a heavy ion collision 
flow as well as it seems to do? A satisfactory matching between the CGC and hy¬ 
drodynamics implies that there should be a certain range of time in which the two 
descriptions predict the same evolution. At the moment, we are not there yet, even 
though considerable progress has been made in the past years. 

After a brief account of why heavy ion collisions are interesting from the point of 
view of QCD (section 2), the rest of this review follows the time-line of a collision. 
We recall the main physical ideas behind the parton model in the section 3, and we 
describe gluon saturation in the section 4. The section 5 is devoted to the Color 
Glass Condensate, the QCD-based effective theory for the saturated regime, and in 
the section 6 we show how to apply it in order to make leading order calculations 
in heavy ion collisions. Next-to-Leading order contributions are considered in the 
section 7, as well as the scale evolution of the gluon distribution in the projectiles 
via the JIMWLK equation. In the section 8, we first introduce the main issues 
and puzzles posed by attempts to match CGC calculations and hydrodynamics. In 
the section 9, we discuss the instabilities that exist in the solutions of the classical 
Yang-Mills equations, and their disastrous consequences for fixed loop-order CGC 
predictions beyond leading order. A resummation that cures these pathologies is 
presented in the section 10, leading to a scheme known as the Classical Statisti¬ 
cal Approximation. We also present alternative derivations of this approximation 
scheme in order to make connections with other approaches. We present in the 
section 11 the results obtained by using this approximation in the context of heavy 
ion collisions, using different types of initial conditions. A discussion of some known 
shortcomings of this approximation is presented in the section 12. 


2. Heavy ion collisions 
2.1. Reminder on QCD 

Although they occupy only a tiny fraction of the volume of atoms, the atomic nuclei 
make up for most of the mass of ordinary matter. The protons and neutrons that are 
contained in nuclei each contain three valence quarks, that give them their quantum 
numbers. However, these valence quarks account only for a small part of the nucleon 
mass. Most of it comes from binding energy, i.e. from the cloud of gluons and virtual 
quark-antiquark pairs that surrounds the valence quarks. This predominance of 
binding energy in the mass of hadrons reflects a crucial property of the force which 
is responsible of the cohesion of hadrons and nuclei: this force becomes strong on 
distance scales comparable to the proton size, around 10~ 15 m. On the other hand, 
the measurements of structure functions in deep inelastic scattering experiments, 
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Fig. 1. Atoms and nuclei. 

first performed at SLAC in the 1960’s, can be understood if one assumes that this 
force becomes weak on distance scales that are much smaller than the proton size. 

The combination of these two properties paved the way to the development 
of Quantum Chromodynamics (QCD), as the microscopic theory that governs the 
interactions between quarks and gluons. On the surface, QCD is a gauge theory that 
resembles very much Quantum Electrodynamics. The matter degrees of freedom 
are spin 1/2 quarks, that interact by the exchange of vector particles, the gluons. 



The quark-gluon coupling in QCD is very similar to the electron-photon coupling 
in QED, except that it has more “structure” since it involves a matrix of the fun¬ 
damental representation of SU(3), t'/. In this object, the index a (running from 1 
to 8, the dimension of the SU(3) algebra) is the color charge of the gluon, and the 
indices i and j (running from 1 to 3, the dimension of the matrices in the funda¬ 
mental representation of SU(3)) are the color charges of the incoming and outgoing 
quarks. The fact that the gluons themselves carry a color charge is the essential 
difference between QCD and QED, since it leads to novel interaction vertices that 
involve only gluons. These new interactions are a requirement of gauge symmetry, 
and can be derived from the following gauge invariant Lagrangian 

£ = - m /)V’/ • ( 2 ) 

/ 

At the classical level, the only free parameters in QCD are the quark masses rrif 
and a coupling constant g. In the quantized theory, the coupling is usually traded 
for a scale a A QCD that emerges from the renormalization of the coupling. This new 

a Note that, in the absence of quarks (or with only massless quarks), QCD is scale invariant at the 
classical level. Loop corrections induce a breaking of this scale invariance, which is the reason for 
the appearance of Aq CD in the quantized theory. 
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scale arises in the running of the coupling constant a s = g 2 /(An). At one loop, this 
is given by 4-8 

t 2nN c (o ^ 

^ ] = (UN c ^2N f )log(E/A QCD ) ’ ( } 

where E is the energy scale, N c the number of colors and Nf the number of quark 
flavors. The main difference compared to QED, due to the self-interactions of the 



Fig. 2. Running coupling in QCD. 

gluons, is the fact that the coupling becomes smaller at short distances as shown in 
the figure 2, a property known as asymptotic freedom. 

A related property of QCD is the long distance behavior of the interaction poten¬ 
tial between a quark and an antiquark. This can be calculated numerically in lattice 
simulations for heavy quarks (that are therefore static). This potential, shown in 
the figure 3, behaves as a standard A/r Coulomb potential at short distance, but 
increases linearly at large distance, in sharp contrast with electromagnetic interac¬ 
tions. This leads to color confinement, that is the fact that free color charges cannot 
exist in Nature. Quarks only appear in color singlet bound states called hadrons, 
made of 3 quarks (baryons) or quark-antiquark pairs (mesons). The spectrum of 
these bound states can in principle be determined from first principles from the 
QCD Lagrangian, and it depends only on the quark masses and on the QCD scale 
A qcd . However, this dependence is non-perturbative and lattice simulations are the 
only way to perform these calculations. Presently, lattice calculations can reproduce 
the spectrum of light hadrons with an accuracy of the order of 5%, as illustrated in 
the figure 4. 

The QCD running coupling shown in the figure 2 can also be viewed with a 
different perspective: it suggests that if one squeezes many hadrons in a small 
volume, then the average inter-quark distance will be small and their interactions 
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Fig. 3. Coulomb potential of a heavy quark and antiquark pair, from lattice QCD. 



Fig. 4. Hadron spectrum from lattice QCD. 

will be weak. In such a situation, the quarks would not be confined into individual 
hadrons, and would instead form a plasma made of deconfined quarks and gluons. 
This idea is substantiated by lattice calculations of the QCD partition function as 
a function of temperature, that indicate a rapid increase of the number of effective 
degrees of freedom at a temperature around 160 MeV b . This suggests that the 
relevant degrees of freedom are no longer the color singlet light hadrons (pions, 
kaons,...) and have been replaced by quarks and gluons (that are more numerous 
because of the uncovered color degree of freedom). 

b This is for QCD with 3 light quark flavors. The transition temperature is higher for pure glue 
QCD. 
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2.2. Heavy ion collisions 

Experimentally, the conditions of such a transition can be realized by colliding 
heavy nuclei at high energy. Such experiments are presently being performed by 
the RHIC (gold nuclei collided at 200 GeV) and by the LHC (lead nuclei collided at 
5.5 TeV). Just after the impact of the two nuclei, the energy density reaches values 
that are more than ten times the normal nuclear matter density, well above the 
energy density at the deconfinement transition inferred from lattice calculations. 
Such a collision, whose total duration is of the order of 10 frn/c, can be divided into 



—> kinetic theory 
—> viscous hydro 
—> classical dynamics 


Fig. 5. Successive stages of a heavy ion collision. 

several stages, as shown in the figure 5. In this figure, we have also indicated what 
kind of tool one may employ for each of these stages. It turns out that macroscopic 
descriptions such as relativistic hydrodynamics are quite successful at describing 
the bulk evolution of the system. Somewhat surprisingly, the matter produced in 
these collisions seems to behave almost like a perfect fluid with close to no viscosity. 
A very small viscosity suggests that this matter is not the siege of strong dissipative 
processes that would rearrange its microscopic degrees of freedom. 

In these lectures, we will be primarily interested in the beginning of the collision, 
up to the point where a hydrodynamical description may become plausible. We 
will adopt a weak coupling perspective 0 , and we will try to follow a heavy ion 
collision in a description which is as close as possible to QCD. Indeed, in collisions 
at very high energy, the initial energy density is so large that the early stages of 
such a collision should be amenable to a weak coupling description, thanks to the 
asymptotic freedom of QCD. Note that a small viscosity, such that could explain 
the success of hydrodynamics, is more naturally obtained in the strong coupling 
limit because the viscosity is inversely proportional to the scattering cross-section 
of the quarks and gluons. However, it is also be possible to get strong interactions 

c In the discussion on gluon saturation, we will see that the emergence of the saturation scale , that 
increases with the collision energy, justifies this assumption. 
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at weak coupling, provided that the occupation number is inversely proportional to 
the coupling g 2 . In this case, the coupling disappears from the scattering rate, and 
the system has many of the features of a strongly coupled system. 

3. Parton model 
3.1. Kinematics 

As discussed earlier, free quarks and gluons do not exist in normal nuclear matter. 
Instead they are confined into color singlet bound states, whose spectrum depends 
non-perturbatively on the parameters of the QCD Lagrangian. The same is true of 
the energy levels of a nucleus: they could in principle be derived from the underlying 
QCD dynamics, but this is even more complicated than in the case of light hadrons 
and at the moment far out of reach of lattice computations. 

Does this mean that we should give up any hope of using QCD to describe 
collisions between such objects? Fortunately, the answer is no, for collisions at 
sufficiently high energy. The kinematics of these collisions is the key to overcome 
this difficulty. Let us consider first a nucleon at low energy (i.e. when the nucleon 
is almost at rest in the observer’s frame), shown in the figure 6. In this cartoon, the 



Fig. 6. Dynamics of the constituents inside a slow nucleon. 


thick lines represent the three valence quarks, and the horizontal axis represents 
time. Only gluon constituents are shown, not the sea quarks. In such a frame, 
the valence quarks orbit with a period comparable to the proton size (they are 
ultrarelativistic). These quarks exchange gluons that provide the binding force, 
which also happens on scales of the order of the proton size. Moreover, the quarks 
and gluons can briefly fluctuate: for instance, a quark can temporarily become a 
quark+gluon state. The lifetime of these virtual states can be anything smaller 
than the proton size d . When studying reactions involving hadrons, one should 
compare the typical timescale of the collision (shown as a blue strip in the figure) 
with the timescales of the internal dynamics of the nucleon. In collisions involving 
low energy hadrons, the hadron has a complicated internal dynamics on timescales 

d But since QCD is a renormalizable theory, the physics of the strong interactions at hadronic 
energy scales does not depend on what happens on much higher energy scales. Therefore, these 
short lived fluctuations have essentially no relevance in hadronic physics. 
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comparable to the duration of the collision, which makes these collisions untractable 
in perturbative QCD. 

Contrast this with what happens in a collision at very high energy. Although 
scattering amplitudes are boost invariant and may be discussed in any frame, it is 
convenient to imagine that we do not change the momentum of one of the hadrons, 
and that all the energy increase is achieved by boosting the second hadron. This is 
illustrated in the figure 7. The blue strip, unchanged compared to the low energy 
case, may be viewed as the size of the first hadron, that we did not boost. All the 
changes are in the internal dynamics of the second hadron, whose timescales are 
now stretched by Lorentz time dilation. The gluon exchanges between the valence 



Fig. 7. Dynamics of the constituents inside a highly boosted nucleon. 

quarks are now happening over timescales that are much longer than the duration 
of the collision, which means that the constituents of the nucleon can be viewed as 
free during the collision. The same happens to the fluctuations of the constituents. 
The lifetime of these virtual states is increased much beyond the collision timescales, 
making these off-shell constituents undistinguishable from on-shell particles®. Since 
there are fluctuations at arbitrary small timescales in a nucleon at rest, increasing 
the energy will uncover more and more of these fluctuations. These simple kinemat- 
ical considerations are the essence of the parton model, 9 that approximates a high 
energy nucleon or nucleus as a collection of quasi-free constituents (called partons), 
whose density grows with energy. 


3.2. Factorization 

From this discussion, it seems that a QCD description of high energy collisions 
between hadrons may be feasible, provided we can provide “snapshots” of their 
partonic content at the time of the collision. What information is necessary in this 
snapshot is not completely obvious at this point, and may vary depending on the 
observable one intends to calculate, but for instance one may think of the following: 

e The concept of on-shell or off-shell particles depends on the duration of the measurement. The 
only way to know that a particle is exactly on-shell is to perform an infinitely long measurement. 
Indeed, an off-shell particle may be viewed as a particle of momentum p whose energy differs from 
the on-shell energy E p (given by the dispersion relation). A measurement that lasts At can only 
resolve energy differences of order 1/At or larger. 
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• flavor and color of each parton, 

• transverse position and longitudinal momentum. 

Of course, these properties of a hadron cannot be known event by event, which 
means that at best a probabilistic description may be achieved, that would allow to 
compute expectation values for event averaged observables. However, the possibility 
of describing hadronic collisions with only a probabilistic partonic description of the 
incoming hadrons is highly non-trivial, because it is an approximation that amounts 
to discarding certain quantum interferences. Without doing any approximation, the 
transition probability from a pair of hadrons hih 2 to some final state X is obtained 
by summing all the relevant reaction channels before squaring the amplitude, 

transition probability _ 1 ^ Amplitudes 2 

from hadrons to X I 2 —^ h\h 2 — > X 

In contrast, the parton model as described above approximates the transition prob¬ 
ability as follows, 

transition probability probability to find I ^ Amplitudes 2 

from hadrons to X ~ pa ^ ns {q, g} in {h 1 ,h 2 } I ^ {q, g} -A X 

ii,a} 

which is clearly not equivalent to the previous formula. This approximation is 
called initial state factorization. Roughly speaking, the physical motivation for 
such a factorization is that the neglected terms are interferences between a hard 
process that occurs on the timescale of the collision and a process internal to one 
of the projectiles, happening on much longer timescales. The vast separation in 
their timescales is what makes the corresponding interference small. At a more 
formal level, this factorization can be established in QCD, with various degrees of 
sophistication 1 depending on the observable. 

3.3. Single parton distributions 

The most developed framework for this type of factorization is the DGLAP formal¬ 
ism, in which one describes the incoming hadrons by single parton distributions. 
These distributions depend on the hadron and on the parton under consideration, 
on the fraction x of longitudinal momentum carried by the parton, and on a momen¬ 
tum scale Q that can be viewed as the inverse of the transverse spatial resolution 
with which the hadron is probed. In the figure 8, the single parton distributions 
of a proton, extracted from deep inelastic scattering data, are shown at a fixed 
resolution scale Q. Although these distributions are non perturbative and cannot 

f The weakest of these factorization theorems are based on leading log factorization, where the 
two formulas are shown to be equivalent for an infinite series of terms of the form (o:. s log(Q)) 71 
(where Q is some hard scale), but not for terms of the form a a (a 3 log(Q)) 71 . Next-to-leading log 
factorization extends the proof of this equivalence to include terms in a s (a s log(Q)) n , and so on. 
All-orders factorization theorems 10 ^ 1 - prove that the two formulas are equivalent up to terms that 
decrease as inverse powers of the hard scale. 
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Fig. 8. Parton distributions of a proton, at the resolution scale Q 2 = 10 GeV 2 . From. 13 


be computed easily from the QCD Lagrangian 8 , QCD predicts how they change 
if one increases the resolution scale Q , via the DGLAP equation. 14 ^ 11 From this 
figure, one sees that the valence quark distribution is predominant at large values of 
the momentum fraction x > 0.1, and is totally negligible at small x. At any value 
x < 0.1, the gluons are the dominant species of partons, and their density increases 
like a power of 1/x when x —> 0. The sea quarks follow the trend set by the gluons, 
but with a suppression factor of order a s since they are produced by the process 
9 -> 99- 

4. Gluon saturation 
4.1. Dense regime of QCD 

Since the DGLAP factorization framework is based solely on the single parton distri¬ 
butions, it is expected to become insufficient at large parton densities. The problem 
that will arise in this regime is illustrated in the figure 9, that shows side to side a 
typical scattering process in the dilute (left) and dense (right) regimes. In the dilute 
situation, the incoming hadrons are “mostly empty”, and hard scatterings are rare 
processes. Moreover, reactions involving more than one parton in each projectile are 
extremely rare (their rate scales as the square of the probability to find a parton). 
But when the parton density is large, processes initiated by multiple partons become 

s Lattice QCD can be used to evaluate the first few Mellin moments of parton distributions, since 
they are given by expectation values of local operators that can be evaluated in the Euclidean 
theory. 
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Fig. 9. Typical partonic processes in a collision between dilute (left) and dense (right) projectiles. 

more likely to happen. Moreover, the partons that lie nearby in phase-space may be 
correlated, and therefore multiparton distributions cannot be inferred simply from 
single parton distributions. A framework that would enable one to calculate these 
processes should provide information about multiparton distributions in hadrons 
and nuclei, and thus should go beyond the DGLAP framework. Moreover, when 
the parton density becomes of the order of the inverse coupling 1 /g 2 , a strongly 
interacting regime -called gluon saturation 1 ^ 3 - is reached, where an infinite series 
of Feynman graphs contribute at each order in g 2 . 

A hint of the fact that the small x saturation regime is qualitatively different 
from the dilute regime appears when plotting the deep inelastic scattering cross 
section slightly differently. This cross-section depends on two independent Lorentz 
invariant quantities, x and the 4-momentum squared Q 2 of the photon exchanged in 
the scattering. However, when plotted against the combination x 0 32 Q 2 1 this data 
appears to line up on a unique curve (see the figure 10). This scaling signals the 



T 


Fig. 10. Geometrical scaling in the DIS cross-section at small x. The horizontal axis represents 
the variable r = x 032 Q 2 . 
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emergence of an x dependent momentum scale, that behaves roughly as Q 2 (x) ~ 
£- 0.32 This scale, known as the saturation momentum, appears as a consequence 
of the nonlinear interactions among the gluons, that become important at high 
density. 

4.2. Saturation condition 

To understand the onset of gluon saturation, it is instructive to go back to the dilute 
regime at large x. In this situation, a hadron appears as a loose collection of a 
few partons. When the hadron is progressively boosted, these partons radiate more 
gluons by bremsstrahlung h , as illustrated in the top panel of the figure 11. As long as 
the density remains low enough, these cascades of gluons develop independently and 
the evolution * 1 of the hadron structure is governed by the linear BFKL equation. 18,19 
Since these additional gluons are contained within the geometrical volume of the 



Fig. 11. Gluon cascades in the small x evolution of a hadron. Top: dilute regime. Bottom: onset 
of the recombination corrections. 

hadron, their density increases rapidly. At some point, their wavefunctions start to 
overlap and their interactions are no longer negligible. Gluons from two different 
cascades can recombine, which tames the growth of the gluon density. Moreover, 
this recombination process makes the x evolution of the gluon distribution non¬ 
linear. 

Before going into a more quantitative description of gluon saturation, it is easy to 
derive a simple criterion for the onset of saturation. Gluon recombination becomes 
likely when the product of the number of gluons per unit area with the cross-section 
for recombining two gluons into one becomes larger than one, 

OsQ^ x A~ 2 / 3 xG {x, Q 2 ) > 1 . (4) 

a gg-*g surface density 

h As discussed before, these gluons are not truly on-shell, but can be viewed as real gluons if the 
lifetime of the quantum fluctuation that gave them birth is longer than the observation time. 

1 Although this terminology is commonly used, it is somewhat of a misnomer, since the hadron does 

not truly “evolve”. It is the observer’s view of the hadron content that changes as the observer’s 

frame is increasingly boosted with respect to the hadron. 
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This condition can be rearranged in order to obtain an inequality on Q : 

Q 2 <Q 2 S = a sxG{x^Q 2 s ) ~ A i/z x - 0 a 

S -v--' 

saturation momentum 


( 5 ) 


This argument justifies the emergence of the saturation momentum, which charac¬ 
terizes the physics of gluon saturation. Saturation is important when the typical 
momentum scales in a process are smaller than Q s . The region where this con¬ 
dition is satisfied is shown in the figure 12. From the more quantitative plot on 
the right panel of this figure, a typical value to keep in mind is that Q 2 is in the 
range 2-4 GeV 2 for nuclei at the energy of the LHC. In the saturation domain, 




Fig. 12. Saturation domain in the x and Q plane. The 3-dimensional plot on the right adds 
information about the A dependence. From. 20 


non-linear gluon interactions are important, which arguably makes the calculation 
of such processes more complicated. However, this also has an unexpected positive 
side: Q s now supersedes all the softer momentum scales in determining the typi¬ 
cal momentum of the relevant partons, and thus also controls the running of the 
coupling. Since Q s increases when x decreases (i.e. when going at higher energy), 
this opens an avenue for an ab initio weak coupling treatment of multiparton inter¬ 
actions (sometimes called the “underlying event” in other contexts) in high energy 
hadronic scatterings. From eq. (5), one sees that the saturation momentum also 
increases with the mass number of nuclei. This implies that, at a given energy, 
saturation effects are stronger in nucleus-nucleus collisions, for large nucleih This 
is important for heavy ion collisions at the RHIC and the LHC, because in these 
collisions the bulk of particle production is controlled by saturation physics. 


] For gold or lead nuclei, the factor A 1 ' 3 that appears in Q 2 is approximately 6. 
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5. Color Glass Condensate 
5.1. Degrees of freedom 

The Color Glass Condensate 14 (CGC) is a QCD based effective theory whose aim is 
to describe quantitatively the gluon saturation regime. The CGC exploits the high 



- .A 1^ 

4 r T ^ n J 


W[p] 


Fig. 13. Degrees of freedom in the CGC effective theory. 

energy kinematics in order to simplify the description of the non-perturbative va¬ 
lence partons. The main idea was already encountered in the qualitative discussion 
of the parton model: a high energy hadronic collision is so brief that the fact that 
the partons are strongly bound by confinement is totally irrelevant. In fact, over 
such short timescales, the internal motion of the partons inside the hadron appears 
completely frozen. Thus, one may view the partons as static in the transverse plane, 
with a large longitudinal momentum. 26-28 For an observer sitting in the center of 
mass frame of the collision, the only information that matters about these partons 
is the color current that they carry along the beam direction. The dominant 
component of this 4-vector is the longitudinal one. In light-cone coordinates, for a 
hadron moving in the +z direction, it reads 

Jai x ) = Pa{x~,x ± )S fl+ , (6) 

where the function p a {x) is the density of color charges of the partons. It does not 
depend on the light-cone “time” x + because of time dilation. Moreover, due to the 
concomitant Lorentz contraction, its x~ dependence is very peaked around x~ = 0. 

It is important to realize that this drastic simplification cannot be used for all 
partons: it is applicable only to those partons whose longitudinal momentum (in 
the observer’s frame) is large enough. The partons that have a rapidity close to the 

k For more detailed reviews of the Color Glass Condensate, one can consult Refs. 21-25 
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observer’s rapidity have comparable transverse and longitudinal momenta, and thus 
cannot be approximated by a longitudinal current. Moreover, for these partons, the 
Lorentz boost factor that slows down their time evolution is not large and one 
cannot neglect their dynamics. Therefore, these slower partons must be treated as 
full fledged quantum fields. 

The situation is summarized in the figure 13: a cutoff y cut must be introduced 
somewhere between the rapidity of the observer and the rapidity of the hadron 
under consideration. The partons close to the observer (mostly gluons, at least at 
leading order in a s ) are described as gauge fields, while those that are close to the 
projectile are approximated as a static longitudinal color current. 

5.2. CGC effective theory 

Thanks to the rapidity separation between the slow and the fast degrees of freedom 
in the CGC, their coupling can be treated as eikonal, i.e. via a term of the form 
J^A^. Therefore, the CGC can be summarized by the following effective action 



(For a collision of two hadrons, the current is the sum of two terms, one for 
each hadron.) The function p a (x~ ,x±) that appears in the current reflects the 
particular arrangement of the fast partons at the time of the collision. It is not a 
quantity that can be predicted event by event, and one can only have a statistical 
knowledge of this object. Therefore, the CGC also introduces a probability distri¬ 
bution W\p\. As we shall see later, all observable quantities must be averaged over 
all the possible configurations of p, according to the distribution W\p\, 



( 8 ) 


In words, one should first calculate the observable for an arbitrary configuration of 
the color charge density p (in a collision of two hadrons, there is a p\ and a P 2 ), 
and then perform a weighted average over all the possible p’s. The justification of 
this procedure will be given in the following two sections. 

6. CGC at Leading Order 
6.1. Power counting 

So far, we have not assumed anything about the magnitude of the color charge 
density p a that describes the fast partons in the CGC effective theory. For the 
CGC to be applicable to the saturated regime, we must allow p a to be as large 
as the inverse coupling 1/g. Indeed, the recombinations due to non-linear gluon 
interactions can stabilize the gluon occupation number at a value of order 1/g 2 , 
where the gluon splittings and the recombinations balance each other. Since the 
occupation number is quadratic in the gauge field, such a value corresponds to 
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p ~ < 7 _1 . As we shall see, such a large value of the source p simplifies the dynamics 
by making it classical at leading order, but complicates things by making an infinite 
set of graphs contribute at each order in g 2 . In order to see this, let us first examine 
the power counting in the CGC effective theory. Consider a generic connected 
graph 1 , as shown in the figure 14. For such a graph, one finds that the order of 



Fig. 14. Generic connected graph in the CGC effective theory. The dots represent the sources p. 
The lines terminated by arrows are the gluons produced in the final state. 


magnitude depends only on the number of produced gluons and the number of 
loops, via the following formula 29,30 

J_ g # produced gluons ^2(# loops) tg\ 

g 2 

The main consequence of assuming that p ~ g~ 1 is that this power counting does 
not depend on the number of sources that are included into the graph. The reason 
for this is that each additional source (of order g~ x ) is attached to the rest of the 
graph by a vertex (of order g ), and therefore does not contribute to the overall 
magnitude of the graph. This also means that, at each order in g 2 , one must sum 
an infinite set of graphs. 

For instance, the inclusive gluon spectrum has the following expansion in powers 
of g 2 


dNi 

d 3 p 


co + ci g 2 + c 2 g 4 d- 


( 10 ) 


where each of the coefficients Co, Ci, • • • is itself an infinite series of terms of the form 
(. 9P ) n , 


ci = ^2ci,n(gpi, 2 ) n - (n) 

n =0 

At this point, we should make an important remark regarding exclusive versus 
inclusive observables. From the above power counting, we see that the average 
number of produced gluons in a high energy nucleus-nucleus collision is of order 


'Typical graphs are not connected: they are made of many disconnected subgraphs. But it is 
sufficient to discuss the properties of one of these subgraphs. 
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1 /g 2 . If we assume for simplicity that the multiplicity distribution is Poissonian 111 , 
the probability to have a given final state (e.g. a final state with a prescribed 
number of gluons) is exponentially suppressed by a factor exp(—This factor 
may be viewed as a Sudakov factor that arises from excluding all the other final 
states. It turns out that these exclusive quantities are very difficult to calculate. In 
particular, the (disconnected) graphs that involve spectator partons contribute to 
exclusive observables, because these spectator partons may end up producing the 
unwanted final states. 

In contrast, many simplifications occur in the calculation of inclusive quantities, 
that involve an average over all possible final states 

(0)= £ V(AA^f)0[f\. (12) 

all final 
states f 

In particular, the high energy factorization results that will be discussed in the next 
section can only be established for these inclusive quantities, and their proof fails if 
one tries to generalize it to exclusive quantities. 

6.2. Calculation of inclusive observables 

The definition of eq. (12) suggests an elementary method for calculating inclusive 
quantities: compute the exclusive probabilities to end up in a given final state /, 
and sum over all possible /’s. However, it turns out that one can calculate them 
in a much more effective way without having to perform explicitly the sum over 
the final states. For this, one should use the Schwinger-Keldysh formalism, 31 ’ 32 in 
which this sum is already “built in”. Any contribution to eq. (12) is the product 
of an amplitude, a complex conjugate amplitude going to the same final state, and 
the observable evaluated on this final state, as illustrated in the figure 15. The 



/ 


Fig. 15. Illustration of the Schwinger-Keldysh formalism. 


diagrammatic rules for the amplitude (right of the dotted line) are the usual time- 
ordered Feynman rules. The propagator is the usual Feynman propagator, which 
for a (massless) scalar particle reads 


G° ++ (P) 


i 

p 2 + ie 


(13) 


m This is not exactly true in the CGC, but this fact does not change the essence of this argument. 
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For the complex conjugate amplitude (left of the dotted line), one needs the complex 
conjugate of the vertices and propagators. The propagator is therefore 

G°__{p) = ^—. (14) 

jr — le 

Across the dotted line, one must use special propagators that represent the on-shell 
particles of the final state /, 

G° + _(p)=2ne(-p°)S(p 2 ) . (15) 

The Schwinger-Keldysh formalism amounts to the following: 

• Draw all the graphs AA —> AA that have a given order in g 2 (the power 
counting is the same as before, with each loop adding one power of g 2 ). 

• Sum over all the possibilities of assigning the labels + and — to the internal 
vertices. 

• Only connected graphs contribute, because when summed over the + and 
— labels, the subgraphs that are not attached to the observable vanish. 


These rules will automatically provide the sum over final states that was included 
in the formula (12). Note that when used in this context, the Schwinger-Keldysh 
formalism is equivalent to Cutkosky’s cutting rules, 33,34 that were developed as 
a tool to compute the imaginary part of transition amplitudes. The superficial 
description of the Schwinger-Keldysh formalism that we have given here can be 
made more rigorous by writing the generating functional for its Green’s functions. 
It can be obtained as follows from the generating functional Z[j ] of time-ordered 
perturbation theory : 


Z\j+,j-\ = exp 


ld 4 xd 4 y G° + _{x,y) 


5j+(x)Sj-{y)_ 


Z[j + ]Z*[j-) 


(16) 


This formula makes more obvious the fact that the Schwinger-Keldysh formalism 
is made of two copies of the ordinary Feynman perturbation theory (one of them 
complex conjugated), “stitched” together by the on-shell propagators G[ j__. 

At this point, we have not really simplified the calculation of eq. (12). It has 
just been rephrased in a more systematic language. The simplifications come from 
noticing the following identities, 


G++ + G — = G + _ + G_+ 

G ++ — G_|_= G_— G_= G r (retarded propagator) . (17) 

When using the Schwinger-Keldysh formalism to calculate inclusive observables at 
leading order, the sum over the + and — indices always generates the combinations 
of propagators that appear in the second of these equations, and these propagators 
therefore become retarded propagators. 
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6.3. Classical equation of motion 


This simplification is particularly dramatic for inclusive observables at leading order. 
The starting point is a double sum over all possible tree diagrams (they all have 
the same order in g 2 when the source is p ~ 3 _1 ) and over all the indices + and 
— of the Scliwinger-Keldysh formalism. Thanks to the previous remark, the second 
sum merely replaces all the propagators by retarded propagators. One is thus left 
with a sum over all the tree diagrams built with retarded propagators, whose first 
few terms would be 



(here for a (f > 3 scalar field theory.) It is then easy to see that this sum is the solution 
of the classical field equations of motion that vanishes when x° —> — oo (this retarded 
boundary condition follows from the fact that we are summing trees that are made 
of retarded propagators 11 ). Although in interacting theories the classical equation of 
motion is a non-linear wave equation, this is a considerable simplification because 
we have now a problem that can be solved numerically. 

The same simplifications work in the case of the CGC: at leading order, it is 
sufficient to solve the classical Yang-Mills equations with null retarded boundary 
conditions 


[D^ F^] = Pl S v+ + p 2 5 V ~ , lim A»(x) = 0 . (18) 

x® — y —OO 


Assuming that we have solved this equation, all the inclusive observables at leading 
order can be expressed in terms of its solution A For instance, the single inclusive 
gluon spectrum is given by 


d/Vi 

dYd 2 p± 


LO 



e i P . {x - y) UxUy A^x)A„(y) , 

A 


(19) 


and the inclusive multigluon spectra simply read 

dN n dN\ 

d 3 P\ ■ ■ ■ d 3 p n lq d 3 pi 

Similarly, the components of the energy-momentum tensor have simple expressions 
in terms of the classical chromo-electric and chromo-magnetic fields E 1 and B l , 

T °o = \ [E 2 + B 2 ] =[ExB] 1 (21) 

fiij 

= — [E 2 + B 2 ] - [E'E 3 + B‘B ] ] . (22) 

n One can see here why it is important to consider inclusive observables for this simplification 
to happen. It is the sum over the final states that leads to the sum over the indices + and 
—. Without this sum, one would be left with time-ordered propagators, which would make the 
boundary conditions of the classical solution untract able. 


dN i 


d 3 Pn 


( 20 ) 
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6.4. Numerical implementation 



Fig. 16. Left: space-time structure of the classical gauge field A 1 ' . Right: 3-dimensional cubic 
lattice. 


Iii order to perform this calculation in practice, one should be aware of the 
following : 


i. High energy collisions are nearly invariant under boosts in the longitudinal 
direction. This invariance has its simplest manifestation if one uses proper¬ 
time (r = V2x+x~ ) and rapidity (rj = | \og{x + / x~)) as the coordinates 
inside the forward light-cone. When written in this system of coordinates, 
the classical Yang-Mills equations do not depend explicitly on rapidity, and 
thus become 1+2 dimensional equations, 
ii. The sources p\ and P 2 have support on the light-cone, where they are 
singular, i.e. proportional respectively to S(x~ ) and <5(a; + ). These sources 
divide the space time in four distinct regions, as shown in the left figure 
16. The gauge field is identically zero in the region 0, and it can be found 
analytically in the regions 1 and 2. 35 In the region 3, the best one can 
do analytically is to obtain the value of the gauge fields and the conjugate 
electrical fields just above the forward light-cone, at a proper time r = 
0+36 . 


= <*\ 


A 0v = 0 


,35 = 0 , o? n = -Ul&Un (n = 1 , 2 ), 
E n 0 = *2 [ a i> a 2] i 


(23) 


where the Wilson lines Ui^(x±) read 


iq f dx (x ,£c i ) 

U 1 (x ± ) = Pe J v i 


( 24 ) 
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Given these remarks, we need only to solve numerically 1+2-dimensional classi¬ 
cal Yang-Mills equations inside the forward light-cone, 37-46 starting with the initial 
conditions (23). Choosing the time variable r determines the Hamiltonian of the 
system, and from there one can determine the Yang-Mills equations in Hamilto¬ 
nian form. In order to handle them numerically, one must discretize space on a 
cubic lattice (see the right figure 16), while time remains a continuously varying 
variable. This is most easily done in the temporal gauge A T = 0 (also called the 




Fig. 17. Left: link variable. Right: elementary plaquette variable. 


Fock-Schwinger gauge in this context). After having adopted this gauge condition, 
the problem has a residual gauge invariance, under any gauge transformation that 
depends only on space. Naive discretizations based on the gauge potentials are 
not very adequate because they lead to violations of this residual gauge invariance. 
Instead, one should adopt Wilson’s formulation, in which the gauge potentials are 
traded in favor of link variables (see the left figure 17), i.e. Wilson lines that span 
one elementary edge of the lattice 

rX+1 

Ui(x) = P expig / ds A z (s) . (25) 

J X 

Under a residual gauge transformation, these links transform as 

Ui(x) -+ Pt(x) Ui(x) fit (a; + i) . (26) 

The electrical fields E l that appear in Hamilton’s equations transform covariantly, 

E\x) -+ Q(x)E i (x)Cl' f (x) , (27) 

and therefore they should live on the nodes of the lattice. In the A T = 0 gauge, the 
Hamiltonian discretized in this fashion reads 
v E'\x)E\x) 

4^ 2 

x\i 

- 4 E 1 - h Re Tr ( u i 0 x ) U i (■* + (* + m ](*) ) • (28) 

Q . o ^ -j 

£c; ij v 

plaquette at the point x in the ij plane 

The only combinations of link variables that enter in this formula are plaquettes 
(i.e. the trace of the product of the four link variables that form an elementary 
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square on the cubic lattice), which are gauge invariant. The Hamilton equations 
that can be derived from this Hamiltonian form a large (but finite) set of ordinary 
differential equations, that can be solved numerically by standard methods such as 
the leapfrog algorithm. 


6.5. Structure of the classical color fields 

At very short times after the collisions (r -C Qj 1 ), the classical chromo-electric 
and chromo-magnetic fields are parallel to the collision axis, 4 ' as illustrated in the 
figure 18. By studying how the expectation value of transverse Wilson loops, 48,49 



Fig. 18. Color flux tubes just after the collision. 


W = exp ig J dx' l A l j , (29) 

depends on the area enclosed by the loop, one can infer the typical transverse size 
of these flux tubes. Indeed, one can roughly view the argument of the exponential 
as the magnetic flux going through the loop. It has been found that W decreases 
roughly as exp (—# x Area) for areas larger than Qf 2 , which indicates that the 
fields are not correlated over transverse distances larger than Qj 1 . Qj 1 is therefore 
the typical radius of these flux tubes. 

From the classical gauge fields, one can compute the spectrum of gluons pro¬ 
duced in a heavy ion collision (see the figure 19). At large transverse momentum, 
the spectrum decreases as fej 4 . Indeed, when k_ j_ Q s , the saturation criterion 
is not satisfied and one should recover the usual perturbative results of the dilute 
regime. In contrast, saturation effects are quite large at small transverse momen¬ 
tum, k± < Q s , where they produce a strong softening of the spectrum, while the 
dilute result would still give a spectrum that grows as Af[ 4 at small k± (because 
there is no dimensionful scale other than k± in the dilute calculation). Let us close 
this section by a remark concerning the energy dependence of the gluon multiplicity, 
obtained by integrating the gluon spectrum over k±. The result is proportional to 
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Fig. 19. Inclusive gluon spectrum. 


the transverse area S± multiplied by Q1, 

iVgluon ~ • (30) 

Q. s 

From this pocket formula, one sees that the energy dependence of the gluon multi¬ 
plicity is directly inherited from that of the saturation momentum, 

iVgluon ~ a ?" 0 ' 3 ~ S° 15 . (31) 


Note that there is no contradiction between the fact that the multiplicity grows 
like a power of the collision energy and the Froissart bound, that tells us that 
hadronic cross-sections cannot grow faster than a ~ log 2 (s). The difference between 
the two is due to the fact that these two objects measure very different things. 
The total cross-section measures the probability that the two projectiles interact. 
As a probability, its growth is constrained by unitarity. In contrast, the gluon 
multiplicity measures the “amount of stuff” which is produced in a collision. It is 
not a probability, and is not bound by the same constraints. More precisely, the 
Froissart bound is related to the growth of the radius of the “black disk” region with 
energy (this radius grows like the log of energy). However, even after a certain region 
has become “black” (and thus its probability of interacting cannot grow anymore, 
because it has reached the unitarity limit), the number of gluons produced in this 
region will continue to grow like a power of energy (because the number of gluons 
per unit area in the incoming projectiles continues to increase like Q 2 ). 
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7. Next to Leading Order 


7.1. Improved power counting 


From the power counting formula derived earlier, we expected for the gluon spec¬ 
trum an expansion of the form 


dN 

d 3 p 


co + ci g +c 2 g + ■ 


(32) 


and we have seen in the previous section how to calculate the term co/g 2 . The 
following terms, Ci,c 2 g 2 , • ■ ■, are given respectively by the sum of 1-loop, 2-loop, 
etc... diagrams. When calculating loops in the CGC framework, we have to recall 
that the degrees of freedom have been divided into sources and fields, separated 
by a cutoff y c u t in rapidity. In the integration over the loop momentum, we must 
use this cutoff in order to prevent the loop momentum to go into the kinematical 
domain which is described in terms of static sources. Failing to do this would lead 
to a double counting of the contribution of the modes that lie in this region. 

In general, loop diagrams will depend on this cutoff: 50 a graph with n loops 
can contain up to n powers of 7/ cu t- Therefore, the coefficients that appear in the 
g 2 expansion of the gluon spectrum can themselves be expanded in powers of the 
cutoff as follows, 


C 1 — Cio T Cn 7/cut 

C 2 = C 20 + C 21 7/cut + C 22 j/cut (33) 

Leading Log terms 

The terms that have the maximal degree in 7/ cut , he. a degree equal to the number 
of loops, are called the leading log termtf’. 

The cutoff 7/cut was introduced by hand when defining the CGC, as a way of 
separating the two kinds of degrees of freedom, and it is therefore not a physical 
parameter. Observables should not depend upon it. As we shall see in this section, 
the leading log cutoff dependence that arises from loop corrections to observables 
can be absorbed into a redefinition of the probability distribution W[p]. This re¬ 
definition turns W[p\ into a cutoff dependent object, but its cutoff dependence is 
universal, which means that the same distributions can be used for all inclusive 
observables. 

Before we continue with a discussion of the leading log terms, let us also mention 
the fact that the next-to-leading log corrections are now known in some cases: in 
the BK equation (a mean field approximation of the JIMWLK equation) 51-54 and 
also for the JIMWLK equation. 55-58 The running coupling corrections have been 
used in some phenomenological studies 59-63 where they appear to be quantitatively 
important. 


°The terminology comes from the fact that y cu t is the logarithm of a cutoff on the longitudinal 
momentum. 
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7.2. NLO result 

Let us give here a sketch of the proof of this factorization. The first step is the 
derivation of an expression for the NLO correction to inclusive observables. It turns 
out that there exists a formal relationship between the LO and NLO contributions, 
valid for any inclusive observable, that reads: 64,65 


O 


NLO 




^(m, v ) + 


a(u)T u 


O, 


(34) 


In this formula, the LO observable 0 LO should be viewed as a functional of the 
initial value of the classical field on some space-like hypersurface (the integrations 
over the variables u and v are on this surface). The operator T u is the generator of 
shifts of this initial condition, in the sense that its exponential translates the initial 
field *4j n it in any quantity that can be expressed in terms of the classical field 

(35) 

The remarkable property of eq. (34) is that the coefficient functions T 2 and a are 
universal: they do not depend on the observable under consideration. Note however 
that although this formula is valid for all inclusive observables, it is not true for 
exclusive observables. 


exp 


o:(M)T t 


7.3. Classical phase-space formulation of quantum mechanics 

In the formula (34), the operator between the square brackets acts only on the initial 
value of the classical fields, while the time evolution from the initial time surface 
to the time where the observable is evaluated remains classical. This is in fact a 
completely general result in quantum mechanics: at the first order in h, the time 
evolution remains classical and h enters only in the initial condition. Let us make 
a digression to justify this important point. This is best viewed if one rewrites the 
evolution equation for the density operator, 

d ^=ih[H,p T \, (36) 

in terms of the Wigner transforms of p T and H , 

s 
2 

I> ■ < 37 > 

(The Wigner transform of H is nothing but the classical Hamiltonian of the system.) 
Note that in these Wigner transforms, the variables x and p are classical phase-space 


W T (x,p) = J ds (x+-\p T \x — 
T-L(x,p) = J ds e l£fr (x + ^\H\x — ■ 
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variables, not operators. It is straightforward to show that eq. (36) is equivalent to 

dW T 2 . / ih, It , , . 

= H(x,p) — sin ( — ( d P d x - d x d P ) J W T (x,p) 

= {n,W T } +0(h 2 ). (38) 

Poisson bracket 

The first line is known as the Moyal equation. It is equivalent to the von Neumann 
equation for p T , except that it is expressed entirely in terms of classical phase-space 
variables. This equation makes the classical limit particularly transparent if one 
expands in powers of H the operator that appears in its right hand side. As one can 
see immediately, its zeroth order in H is the usual Poisson bracket. This means that 
at the order h °, one recovers classical Hamiltonian mechanics. A remarkable feature 
of the Moyal equation is that it has no term of order H 1 . This means that at NLO, 
the time evolution of the system remains purely classical. The only quantum effects 
at order h 1 come via the initial condition, through the fact that the support of the 
Wigner distribution of a quantum state must have an extension of at least K. The 
first correction to the Moyal equation arises at the order fi 2 , i.e. at NNLO. At this 
order, one gets quantum corrections both in the initial condition and in the time 
evolution itself. This discussion also indicates that a formula such as eq. (34), where 
only the initial conditions are altered, presumably does not exist beyond NLO. 


7 . 4 . Cutoff dependence 

Eq. (34) is very useful in order to extract the cutoff dependence in inclusive observ¬ 
ables at NLO, because this cutoff dependence is already present in the operator that 
acts on O lo . If we keep only the terms that are linear in the cutoff, we have 64-66 

\ J r 2 (u, v) + J ol(u) T„ = y+ ut Hi + y~ ut H 2 ■ (39) 


In this equation, and y~ ut are the cutoffs corresponding to the right and left 
moving nucleus respectively, and Hi .2 are operators known as the JIMWLK Hamil¬ 
tonians 9 of the two nuclei, 67-75 




S ^ , 6 

TXab{X±,y±)~ 


where 


Xab(x±,y±) = 


x± ,y ± Spa(£±)' aDK S P b{y±) ’ 

(x ± - z ±) • (y L - Zj) 


(40) 


47T 3 


/ 


d z± 


(x± - z±) 2 (y± - ziff 1 


x \l-U\xffU{zi_))(l-U\zffU{yff 


ab 


(41) 


In this equation, U is a Wilson line in the adjoint representation, constructed from 
the gauge field A + such that V^A + = —p. 

p If one expands this Hamiltonian at small p, one can recover the BFKL equation , 18 d9 
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7.5. Factorization 

The formula (39) has two important consequences: 

i. This formula is the sum of two terms corresponding to the two nuclei, but 
there is no cutoff dependent term mixing the two nuclei. This means that 
the cutoff dependent terms are intrinsic properties of the nuclei prior to 
their collision, and this is the reason why it is possible to eliminate them 
by redefining the distributions of the sources p±^ of the two projectiles, 
ii. Since the operator in the square brackets in eq. (34) is the same for all 
inclusive observables, the cutoff dependence is equally universal. This is 
the reason why it will be possible to define cutoff dependent distributions 
W[p\^] such that they cancel the cutoff dependence of all observables. 

The property i, about the absence of mixing of the cutoff dependence between the 
two nuclei, can be understood simply in terms of causality. This is illustrated in the 
figure 20. Indeed, the cutoff dependence arises from the phase-space integration of 


Tcoll E 



Fig. 20. Causality argument for factorization. 


the soft gluons emitted by bremsstrahlung. Since they are soft, the formation time 
of these gluons is large: they cannot be emitted during the very brief duration of 
the collision, so they have to be emitted before the collision. Because the separation 
between the two nuclei is space-like until the collision, causality forbids any cutoff 
dependent term that would mix the two nuclei. The property ii also follows the 
same reasoning: gluon emissions that happen before the collision should be the 
same for all observables measured after the collision. 

If we compute observables for fixed configurations p\ and p 2 of the color charge 
densities in the two projectiles, there is no way to get rid of the cutoff dependence. 
The only way to remove it is to integrate over all the possible configurations of p\ , 2 - 
The main ingredient in this manipulation is the fact that the JIMWLK Hamiltonian 
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7~L is a self-adjoint operator: 

J[Dp] w (HO) = J[Dp] (HW)O . (42) 

This property can be used to transfer the action of H from the observable onto the 
distribution W\p\. From eq. (39), one can see that p-averaged quantities such as 


dNi 

(Dp Leading Log 


[D Pl Dp 2 \ Wx[ Pl ] W 2 \p 2 \ 


dNi 



fixed pi ,2 


(43) 


are independent of the cutoff, provided that the distributions W[p\ themselves de¬ 
pend on the cutoff according to the JIMWLK equation 


dW 

dy 


= nw. 


(44) 


From eqs. (34) and (39), it is furthermore obvious that the same factorization for¬ 
mula (with the same W’s) applies to any inclusive observable. 

In eq. (43), it is the evolution with rapidity of the distributions W that gives 
the gluon spectrum its rapidity dependence. Indeed, the gluon spectrum for fixed 
Pi 2 that enters in the integrand is still independent of rapidity. From the JIMWLK 
equation, one sees that the distributions W evolve significantly for changes of the 
rapidity of the order of Ay ~ aj 1 . This factorization result thus predicts that the 
gluon spectrum is rather flat in rapidity at weak coupling. 


7.6. Ridge correlations 


Eq. (43) can be generalized to the inclusive multigluon spectrum. Recalling also 
eq. (20), we obtain the following factorization formula 


dN n 


d 3 p 1 • • • d 3 p n Leading Log 

= J [D Pl Dp 2 ] Wi [ Pl ] W 2 [p 2 \ 


dN i dN i 


d 3 pi d 3 p n 


(45) 


This equation, valid at leading log accuracy, shows that at this order the corre¬ 
lations between the produced gluons only originate from correlations of the p’s of 
the incoming projectiles, since in the integrand the n gluons still appear completely 
factorized. Since the relevant rapidity interval for a significant JIMWLK evolution 
is Ay ~ aj 1 , this is also the typical rapidity distance over which the produced 
gluons will be correlated. 

This is the basis of an interpretation of the peculiar shape of the 2-hadron 
correlations observed in heavy ion collisions. This correlation function is represented 
in the figure 21, as a function of the relative azimuthal angle and relative rapidity. 
As one can see, the correlation is very narrow in azimuthal angle, and very elongated 
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Au+Au central 
3<p t trj 9<4 GeV/c 


o 

■°- 5 [^\ 


Fig. 21. Two hadron correlation measured in heavy ion collisions.' 6 


in rapidity q . Because of causality, the existence of a correlation between particles 
that are widely separated in rapidity must originate from phenomena that happened 
very shortly after the collision. This is explained in the figure 22. Let us consider 



detection (~1 m/cj 

freeze out (-10 fm/c) 
latest correlation 


Fig. 22. Origin of the rapidity correlations. 

the time evolution of two particles A and B in reverse, starting from their last 
interaction on the freeze-out surface. Obviously, by causality, they must come from 
the light-cones represented respectively in red and green in the figure. A correlation 
is an event that had an influence on both the particle A and the particle B. It must 
therefore have happened in the overlap between these two light-cones, that we have 
represented in blue. One sees clearly that there is a maximal time at which this 
correlation could possibly have been created. From the time of the freeze-out and 
the rapidity separation of the two particles, it is easy to determine this upper bound 
of the time, 

^correlation A 7"f re eze out C ^ ^ ^ ■ (46) 


q The peak in the middle is a jet-like correlation, due to quasi-collinear splittings in the final state. 
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This bound decreases very rapidly as one increases the rapidity separation A y. In 
heavy ion collisions, the order of magnitude of the freeze-out time is 10 fm/c. For 
instance, a correlation between particles separated in rapidity by A y = 6 must be 
produced before the time 0.5 fm/c, which is well within the regime where the CGC 
is still applicable. 

It is in fact easy to understand qualitatively the main features of the observed 
correlation from the structure of the classical color fields produced at early times 
in heavy ion collisions.' 7,78 As we have said before, these fields are organized in 
flux tubes that have a typical transverse size of Qj 1 and that remain coherent 
over rapidity intervals of order aj 1 . Two gluons emitted from the same tube are 



Fig. 23. Rapidity correlations from color flux tubes. 

correlated if they are produced with a rapidity separation Ay < aj 1 , but are not 
correlated if they come from two different flux tubes r . From the size of the flux 
tubes, we conclude that the probability that two particles are correlated scales as 
( RQ S )~ 2 where R is the transverse radius of the collision zone. 

This explains the existence of a long range correlation in rapidity between pairs 
of particles, but not why their correlation is peaked in azimuthal angle. The 2- 
gluon correlation one gets from the CGC color fields is nearly independent of the 
azimuthal angle, because on average these gluons can be emitted in any transverse 
direction. However, one should keep in mind that the above causality argument 
applies only to the correlation in rapidity, not to the correlation in azimuthal angle 
that can be produced much later. These azimuthal correlations are generated by the 
radial flow' 9,80 that develops subsequently and expels radially the matter produced 
in the collision. Simple relativistic kinematics indeed shows that if one boosts a 2- 
particle spectrum independent of azimuth, it becomes peaked around A0 = 0 (the 
prominence of the peak increases with the velocity of the boost).'' In the figure 
24, we show a comparison of the strength of the azimuthal correlation in data and 
in a very basic radial boost model where a unique radial boost velocity is applied 
to a flat spectrum (the radial velocity is estimated from the slope of p± spectra at 

r The fact that the chromo-electric and chromo-magnetic fields are purely longitudinal at early time 
does not seem to play any role in this argument. The important properties are their coherence 
length in rapidity and in the transverse plane. 
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Fig. 24. Comparison of the strength of the azimuthal correlation in data and in a simple radial 
boost model 81 (see also 82 ). 

small momentum). The centrality dependence, which in this model comes from the 
increase of the radial velocity with centrality, is in fair qualitative agreement with 
the measurement. 

8. From the Color Glass Condensate to hydrodynamics 
8.1. Requirements for hydrodynamics 

The CGC provides a self-contained QCD based framework for describing heavy 
ion collisions from first principles. It also provides tools for calculating inclusive 
observables at leading log accuracy, i.e. leading order plus a resummation of all 
the leading log contributions coming from higher loop diagrams. However, there is 
some physics that plays an important role in heavy ion collisions but is not easily 
captured in the CGC. The fact that the produced gluons and quarks will eventually 
hadronize when their energy density falls below the QCD critical energy density is 
obviously not present in the CGC (at least, at any fixed loop order). 

For this reason, CGC calculations can only describe the early stages of the 
fireball expansion, and should be later on matched onto another description such 
as relativistic hydrodynamics. 83-90 Since the components of the energy momentum 
tensor can be computed in the CGC framework, one could use it as initial data for 
the hydrodynamical evolution. Firstly, for such a matching to be possible, the CGC 
must bring the system to a state that hydrodynamics can handle. This means that 
the transverse and longitudinal pressure should not be too different (in particular, 
the longitudinal pressure should not be negative), and that the viscous effects (e.g. 
the ratio y/s of the shear viscosity to the entropy density) should be small. 

Moreover, when performing such a matching between two descriptions that use 
different degrees of freedom, one should be careful to check that the two descriptions 
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Fig. 25. Smooth matching between CGC and hydrodynamics. 


are compatible in a certain time window. In other words, there should be some range 
where the two models predict equivalent results. If this is the case, the precise time 
To at which the matching is realized is not important, and it can be varied in this 
range without any incidence on the final result s . 

The typical behavior of the ratio ?;/s in a gauge theory in equilibrium is shown 
in the figure 26. When the coupling is small, this ratio can be calculated in a weak 


r| Is 



Fig. 26. 


Shear viscosity to entropy ratio in a gauge theory as a function of the coupling. 


s This is very similar to the factorization of the source distribution W\p\ in the calculation of 
inclusive observables. The independence with respect to the cutoff y cu t is possible because the 
static sources describe the same physics as the gauge fields in a certain range of longitudinal 
momentum. 
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coupling expansion. For QCD, it is given by 


r) ^ 5.1 

s < 7 4 In 



(47) 


This formula shows that 77 /s is large at weak coupling. At large coupling, this quan¬ 
tity cannot be calculated in QCD, but there is an exact result for a supersymmetric 
cousin of QCD, the J\f = 4 SUSY Yang-Mills theory*: 77 /s = l/( 47 r ). 91,92 From this 
plot, it seems that only strongly coupled systems can have a small 77 /s ratio. There 
is however another possibility to evade this conclusion. Firstly, one should recall 
the kinetic interpretation of the ratio 77 /s, 

rj mean free path 

s de Broglie wavelength 

In a system where the degrees of freedom have a typical momentum Q , the de 
Broglie wavelength is of order Q -1 , while the inverse mean free path is given by 


[ fk (1 + fk) 
Jk 


(mean free path) 1 ~ g 4 Q 2 x 


(49) 



cross section density 


Bose 

enhancement 


The factor 1 + fk under the integral is not needed when discussing dilute plasmas, 
but is important if the occupation number is large. In the CGC, just after the 
collision of two heavy ions, one has /*. ~ g~ 2 up to k ~ Q. Therefore, in such a 
system one has 77 /s ~ g°, which is much smaller than the perturbative result in 
equilibrium. It seems therefore plausible that the strong color fields produced in 
heavy ion collisions may flow, not because they are strongly coupled but because 
they are highly occupied. 

8.2. Expansion and free streaming 

The other main feature of heavy ion collisions is the very rapid expansion of the 
system in the longitudinal direction, which causes a redshifting of the longitudinal 
momenta. As illustrated in the figure 27, the pressure tensor tends to become 
anisotropic because of this, unless the interactions are strong enough to overcome the 
expansion. The figure shows what the expansion does on non-interacting particles. 
Starting from a nearly isotropic distribution of the velocities at the time t\, the 
expansion will “filter out” the particles so that at the time T 2 only particles with 
the momentum rapidity y ss rj exist at the space-time rapidity 77 . This means that, 
in the local comoving frame, the longitudinal pressure is much smaller than the 
transverse pressure. A large pressure anisotropy is a problem for hydrodynamics. 
Indeed, the difference between the longitudinal and transverse pressures goes into 

‘This gauge theory is conformal, and is dual to a string theory in an AdSs X S5 background. The 
large coupling limit of the gauge theory corresponds to the weak coupling limit of the string theory, 
in which gravity becomes classical and reduces to Einstein’s equations. 
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Fig. 27. Role of the expansion in decreasing the longitudinal pressure. 


the viscous terms 11 , and large viscous corrections are a sign that hydrodynamics 
may incorrectly reproduce the underlying dynamics. 

Let us now describe the CGC prediction for the energy momentum tensor, start¬ 
ing with the LO calculation. Immediately after the collision, at r = 0 + , it is known 
analytically that the chromo-electric and chromo-magnetic fields are both parallel 
to the collision axis. This peculiar structure of the color fields implies that the 
energy momentum tensor is diagonal, = diag(e, P T ,P T ,P L ), with w P T = e and 
P L = — e. At later times, the energy-momentum tensor must be determined numer¬ 
ically by solving the classical Yang-Mills equations, and by computing T^ from the 
classical gauge fields according to eqs. (22). The results of this calculation 47,100 are 
shown in the figure 28. After starting at —1, the ratio P L /e increases and becomes 
mostly positive at a time of order Q s t ~ 1. However, this calculation shows that 
the longitudinal pressure remains at all times much smaller than the transverse one. 
Thus, the CGC at leading order leads to a situation which is similar to free stream¬ 
ing, where e,P T ~ r^ 1 and P L -C . This leads to an unsatisfactory matching 
between the CGC at LO and hydrodynamics 51 , as illustrated in the figure 29. In¬ 
deed, in hydrodynamics the ratio P L /P T increases to approach 1 while it remains 
near zero in the CGC at LO. 


U There are now attempts to view hydrodynamics as an expansion around a non-isotropic back¬ 
ground. In this formulation, this may be less of a problem. 93-99 

v This discussion also applies to the LO result improved by the resummation of the leading log 
corrections. This does not change the conclusion of this paragraph since this resummation is 
totally absorbed into the rapidity evolution of the distributions W[p]. 

w The negative longitudinal pressure in a longitudinal flux tube is the analogue of a string tension. 
x In principle, this matching requires several steps: (1) compute T^ u {x) from CGC, (2) find its 
time-like eigenvector such that u^T^ Ll/ (x) = eu 1 ' (this defines the local rest frame, and the local 
energy density), (3) compute the pressure from some equation of state P = /(e), (4) compute 
the viscous stress tensor as the difference between the full and the ideal part (obtained from 
e, P and u ^). In many calculations using “CGC initial conditions”, a simplified procedure is often 
employed, where one assumes that u^ = (1,0): (1) compute e = T°°, (2) define P = /(e), (3) 
neglect the viscous stress tensor. 
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Fig. 28. Transverse and longitudinal pressure to energy density ratios, in the CGC at leading 
order. 



Fig. 29. Matching between hydrodynamics and the CGC at LO. 


8.3. A simpler model study case: thermalization in a box 

Before presenting recent works on isotropization and thermalization in heavy ion 
collisions, it is interesting to pause a brief moment on the same question for a 
hypothetical system of gluons enclosed in a fixed volume. In this case, the ultimate 
outcome is completely clear from the start: the system will eventually thermalize, 
and since its total energy is conserved one can predict from the start what will be 
the equilibrium temperature. The only issues are the timescale of the thermalization 
process, its possible dependence on the nature of the initial condition, and the shape 
of the gluon distribution at intermediate times before thermalization is complete. 
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A particularly interesting class of initial conditions are the so called “overoccu¬ 
pied” initial distributions, because they are also realized in heavy ion collisions in 
the CGC framework. The typical CGC-like initial condition has modes occupied up 
to the saturation momentum Q s , with a large occupation number of order g~ 2 . It 
is called overoccupied because it contains too many particles for its energy density: 

Q 3 s Qt 


- 3/4 


r 


r 


9 


- 1/2 


» 1 


(50) 


while the dimensionless combination ne~ 3 / 4 should be of order 1 in a system with 
the same energy density in equilibrium. With such an initial condition, the mean 
free path in the system is parametrically shorter than the thermalization time, 
which leaves ample time for the system to forget its initial conditions long before 
reaching thermal equilibrium. With this type of initial condition, the hard scale of 
the system grows with time according to the following law 101-103 


A h ~ Qs ( t/tof 17 > 


(51) 


and the distribution scales as 


f(t,p) ~ (toA) 4/7 f(p/A H ) . (52) 


The thermalization time can be estimated as the time at which the hard scale 
reaches the expected equilibrium temperature (i.e. the fourth root of the initial 
energy density), T eq ~ <3 s g -1 / 2 . This gives a thermalization time that has the 
following parametric dependence on the coupling, 

Qste q ~ g~ 7/2 . (53) 


In a recent work 104 using the effective kinetic theory y of ref., 105 this estimate was 
made more quantitative and the numerical prefactors computed to give 


72 1 

:eq ~ l + 0.121og(A- 1 ) ’ 


(54) 


where A = g 2 N c . 

This type of overoccupied initial condition has also led to speculations about 
the possibility of forming a Bose-Einstein condensate (BEC) in such systems, that 
would accommodate the excess of particles. 102 Such a condensate can only be tran¬ 
sient, because the number of particles is not conserved in a relativistic quantum 
field theory. Therefore, the true equilibrium state cannot have a BEC. Whether 
such a condensate can form as a transient phenomenon depends on the magnitude 
of the initial overoccupation and the rate of the number changing processes. In a 
scalar </> 4 theory, where the rate of the inelastic processes is quite small compared to 
the elastic one, such a condensate has been observed in a number of numerical sim¬ 
ulations. 106,10 ' The situation is much less clear in QCD, because the inelastic rate 

y This effective theory includes 2—^2 processes dressed by in-medium masses, and effective 1 —> 2 
and 2 —> 1 processes due to the quasi-collinear splitting of hard gluons by bremsstrahlung. For 
the latter, one must resum multiple scatterings, due to the Landau Pomeranchuk Migdal effect. 
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is enhanced by soft and collinear divergences. Unsurprisingly, kinetic theory cal¬ 
culations (based on a small scattering angle approximation of the matrix element) 
neglecting the inelastic processes do observe the formation of a BEC. 108-110 Exten¬ 
sions of these QCD kinetic theory computations to include inelastic processes 111 
suggest that a BEC may still appear after including number changing processes, 
while other computations do not see any evidence for it. 103,104 

Note that for initially underoccupied 2 initial conditions, the evolution towards 
equilibrium is quite different. In particular, it is a lot less universal because the 
mean free path associated with the initial distribution is comparable to or longer 
than its value at equilibrium. In this situation, the initial hard particles first create 
a bath of soft particles by bremsstrahlung radiation, which thermalizes quickly. 
Then the remaining hard particles lose their energy by successive 1 —> 2 splittings 
induced by collisions on the soft background. 

9. Weibel instability and resummation 
9.1. Unstable classical solutions 

There is however a good reason to explore the CGC beyond leading order: namely, 
the fact that the boost invariant solutions of the classical Yang-Mills equations are 
unstable. 100,114-122 When their initial condition is modified by a small but rapidity 
dependent perturbation, the solution diverges from the unperturbed one. This is 
illustrated in the figure 30, that shows a component of the energy momentum tensor 
that should be very small at all times if the perturbation was stable. Instead, it 
growths like an exponential exp y/Jrf (the characteristic growth rate /r is of the 
order of the saturation momentum Q s ). These unstable modes in the classical 
Yang-Mills equations are closely related to the Weibel instability that occurs in 
anisotropic plasmas in QED and in QCD. 101,113 ’ 123 " 142 More details on how these 
instabilities of the classical solutions develop can also be found in re fs. 100,120,121,143 

As we shall see in the next subsection, these instabilities are disastrous for the 
power counting that we have exposed earlier, where one keeps track only of the 
powers of g 2 . Indeed, terms that have a higher order in g 2 because they arise 
at a higher loop order may in fact contain time dependent factors that increase 
exponentially with time. These secular divergences mean that fixed loop order 
calculations are most likely unreliable (we will show an example of this in the next 
subsection), and that the power counting should be revisited and improved in order 
to capture the most important among those terms. 

On the other hand, the existence of instabilities in the classical solutions of the 

z In the bottom up scenario of ref. 112 (see also ref. 113 ), it has been argued that in heavy ion 
collisions the expansion may turn a CGC-like initial condition into an underoccupied distribution 
before full thermalization is reached. 

a The fact that y/r appears here instead of r itself is due to the longitudinal expansion of the 
system. Because of the expansion, the equation that drives the growth of the perturbations is a 
Bessel equation instead of an ordinary wave equation. 











September 7, 2015 1:7 


World Scientific Review Volume - 9.75in x 6.5in 


QGP.5 page 38 


38 


F. Gelis 


0.0001 



le-05 


le-13, 


0 500 1000 1500 2000 2500 3000 3500 


2 


g 


Fig. 30. Growth of unstable modes in classical Yang-Mills dynamics. 


Yang-Mills equations may play a very important role in the statistical equilibration 
of the system. If we consider this system from the point of view of its classical phase- 
space, its initial condition at r = 0 + occupies a very compact region in phase-space. 
Indeed, at LO, the support of its Wigner distribution is a single point, since the LO 
is purely classical. At NLO, this support is broadened into a region of extension h. If 
the dynamics was completely stable (as in an integrable system for instance, where 
the number of conservation laws equals the number of degrees of freedom), the size 
of this support would remain roughly constant in time, and at any subsequent time 
the system would still be far from statistical equilibrium. In contrast, an unstable 
dynamics will map this initially compact support into the extended portion of the 
phase-space allowed by the few conservation laws that remain valid. Measurements 
performed after such an evolution should bring results that are more in line with 
microcanonical equilibrium. 

9.2. Pathologies at fixed loop order 

These instabilities force us to reconsider the power counting that was the basis for 
organizing the expansion in powers of g 2 of inclusive observables. As we have seen 
before, the one-loop corrections -that are formally of relative order g 2 - contain 
leading log terms proportional to the cutoff y c u t, that can be absorbed into a redef¬ 
inition of the distributions W[p]. Because of the instabilities, the 1-loop correction 
also contain some terms that grow exponentially in time. The best place to see this 
is via the following expression for the 2-point function T 2 (u,v) that enters in the 
formula (34), 



( 55 ) 
















September 7, 2015 1:7 


World Scientific Review Volume - 9.75in x 6.5in 


QGP.5 page 39 


39 


(to is the initial time surface on which eq. (34) is expressed) where the functions 
afc(r, x) are small perturbations around the classical field encountered at leading 
order (for Yang-Mills theory, these functions would also carry color, spin and Lorentz 
indices not written explicitly in eq. (55)). They obey the equation 


VpVW-V'jr + igFn 


at = 0 


lim at(x) = e^(fc) e 


ik-i 


(56) 


in which the covariant derivatives contain the LO classical field. The plane wave 
initial condition used in eq. (56) for these small perturbations can be understood 
from the explicit derivation of the formula (34). 

The existence of instabilities means that, for some range of k, the mode func¬ 
tions afc grow exponentially with time. The naive power counting derived so far 
implicitly assumed that the function T 2 is of order g °, but we never looked at its 
time dependence. From eq. (55), it is now obvious that it will become very large 
if some of the a^s are unstable. Eventually, these exponential factors in T 2 will 
compensate the g 2 that comes from the loop, and these terms will be as large as 
the LO terms. This statement is illustrated in the figure 31, where we compare the 




Fig. 31. Effect of parametric instabilities on the perturbative expansion of the energy-momentum 
tensor in a scalar (/> 4 theory. 

energy density and pressure at LO and LO+NLO in a </> 4 scalar field theory, 144 for 
which the solutions of the classical field equations of motion are also unstable 13 . This 
computation shows clearly that the fixed order LO+NLO result cannot be trusted 
after some time, because it becomes much larger than the LO result (note that this 
happens only for the pressure, because the energy density is protected from this 
exponential growth by energy conservation). 

Although a similar NLO calculation has not been done in the CGC framework, 
one can guess what would happen. The Weibel instabilities would produce an 

b The instability in the </> 4 theory is of a totally different nature, since it is caused by paramet¬ 
ric resonance. 145 However, the detailed mechanism of the instability is not important in this 
discussion. 
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unbound growth of the longitudinal pressure, and this time the ratio P L /P T would 
be driven to arbitrarily large values exceeding unity. Attempting to match such a 



Fig. 32. Matching between hydrodynamics and the CGC at NLO. 


NLO CGC initial condition to hydrodynamics would not be better than at LO, as 
illustrated in the figure 32. 

9.3. Resummation of the leading terms 

The same problem would in fact occur at any fixed loop order. The only way to 
improve the situation is to examine each loop order and to keep the most important 
terms at each order. For this, we need first to modify the power counting rules that 
we have established earlier, in order to keep track of the unstable modes. Let us 
first examine the graph that contributes at 1-loop, represented in the top-left corner 
of the figure 33, in conjunction with the formula (34). In this formula, each of the 
derivatives with respect to the initial classical field creates a perturbation to this 
classical field, that we have indicated by green propagators in the figure 33 (in the 
top-left graph, only the term proportional to r 2 , that has second derivatives with 
respect to the initial fields, has been represented - the term with only one derivative 
has a slower growth with time). The “standard” power counting 13 would assign a 
factor 1 to r 2 and a factor g to each of the derivatives T u (represented by a blue 
dot in the graphs). Thus the NLO correction to the energy-momentum tensor is 
expected to be of order g°, while the LO is of order g~ 2 . From this diagrammatic 
representation, it also easy to count the number of perturbations of the classical 
field. Each of them will develop into a factor of order exp (g is of order Q s ). 

c For the term in c*T in eq. (34), the standard counting has a factor g from ex. and a factor g from 
the operator T. However, we will not need to consider this term further since it has a subleading 
growth in time. 
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Fig. 33. Improved power counting taking into account the growth of the unstable modes. 

Thus, the expansion of T 7 *" is more accurately written' 1 as : 

T» v = c 0 g~ 2 + ci g° e 2 ^ + • • • , (57) 

where the coefficients Co, ci do not grow exponentially with time. From this pocket 
formula, one can deduce at which time the naive loop expansion breaks down. This 
is the time when the one-loop result becomes as large as the leading order, i.e. 

Tmax ~ /i -1 log 2 (l/0 2 ) . (58) 

Up to a logarithmic factor, this time is of the order of the inverse saturation mo¬ 
mentum. 

At two loops, the naive power counting tells us that we should get terms of 
order g 2 (i.e. g A relative to the leading order). However, not all the terms have 
the same growth in time. In the figure 33, we have represented two types of 2-loop 
contributions, in order to illustrate these differences. In the top-right graph, the two 
loops are the seed of four perturbations to the classical field, while in the bottom 
graph, only three of these perturbations are created. The latter term will therefore 
have a subleading behavior in time. Moreover, one sees that the distinguishing 
feature of the top-right graph is that it can be generated by acting twice with the 

d This formula indicates the worst possible behavior. It is possible that some components of the 
energy-momentum tensor will not be affected by the instability, as was the case in the scalar field 
theory considered in the figure 31. 
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quadratic part of the operator that appears in eq. (34). This reasoning extends to 
all orders. At the n-tli loop order, the maximal number of field perturbations that 
can be seeded on the light-cone is 2 n, and the corresponding graphs are generated 
by acting n times with this quadratic operator. The sum of all these leading terms 
can be obtained by 


= exn 

resummed * 


^ j r 2 (tt,v)T u T l 


i LO 


(59) 


Note that the Taylor coefficients of the exponential correspond precisely to the 
symmetry factors of graphs such as the top-right diagram of the figure 33 (the 1/2! 
of the second Taylor coefficient gives the symmetry factor that corresponds to the 
freedom of swapping the two r 2 ’s that are hanging below the light cone). 

Moreover, if the 2-point function T 2 used in eq. (59) is precisely the one that 
enters in the NLO result (34), then one has 


rpllf _ rpnv . ('gQ') 

^resummed ~lo ' nlo ' 

In other words, this resummation contains the exact LO and NLO contributions, 
and a subset of all the higher loop contributions. It is important to keep in mind 
that, starting at the 2-loop order, it is an approximation which is not equivalent 
to the complete underlying theory. This will have important consequences that we 
will discuss later. 


10. Classical statistical approximation (CSA) 


10.1. Reformulation as a Gaussian average 


At this point, the resummation performed via the eq. (59) is quite formal. Three 
questions must be addressed: (1) can this formula in terms of functional derivatives 
be evaluated in a practical way? (2) does eq. (59) lead to results whose time 
dependence is bounded? (3) when doing this, are we introducing other pathologies? 

In order to answer the first question, one should recall the following identity: 

+oo 2 

2 t P ~ Z / 2a 

e? 9 * f(x) = dz — . f(x + z) . (61) 

J \f2-Ka 


This formula can be established e.g. by applying a Fourier transform to both sides. 
Although we cast it here in a space of functions of a single variable, this formula 
can be generalized to operators that are Gaussian in derivatives over a functional 
space. It enables us to rewrite eq. (59) as® 


T t—^ = J i Da ] ex P - l J a (“) r 2 \u lV )a(v) 


TglAinit + a] ■ (62) 


e Such a Gaussian averaging procedure has also been reached in other approaches. 146 150 
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This resummation procedure, where one averages classical trajectories over an en¬ 
semble of initial conditions, is known as the Classical Statistical Approximation 
(CSA). 

From this equation, one can easily see that the problem of the unbounded growth 
of the fluctuations has been cured. Indeed, this resummation has promoted the 
linearized perturbations 1 that appear in the NLO contribution into an integral part 
of the non-linear classical field (the initial condition of the classical field is modified 
by the perturbation, but its evolution remains fully non-linear). In any theory where 
the potential prevents the fields from running away to infinity, this guarantees that 
the resummed quantity will not diverge in time. 

10.2. Practical implementation 

In the form of eq (62), the procedure for evaluating the resummed energy- 
momentum tensor 8 is quite straightforward: 

1. Determine the 2-point function T 2 (u,v) that defines the Gaussian fluctu¬ 
ations, for the initial time Q s tq of interest. This is an initial value prob¬ 
lem, whose outcome is uniquely determined by the state of the system at 
x° = —oo, and depends on the history of the system from x° = — oo to 
r = To- This problem is solvable analytically as long as the fluctuations 
remain weak, a M <C Q s /g • If they grow larger, the fluctuations start to 
interact non-linearly, and their spectrum becomes non-Gaussian. To avoid 
this, the initial time should be chosen such that Q s to <C 1. 

2. Solve the classical Yang-Mills equations from To to t/. In high energy 
collisions, the problem as a whole is boost invariant, but individual field 
configurations are now rapidity dependent because of the fluctuating part 
superimposed to their initial condition. Therefore, unlike in the CGC at 
LO, the classical Yang-Mills equations must now be solved in 3+1 dimen¬ 
sions. 

3. Do a Monte-Carlo sampling of the fluctuating initial conditions. 

The setup for doing this is the same as the one described when discussing the 
CGC at LO, except for the fact that one must keep rapidity explicitly. One must 
discretize the classical Yang-Mills equations in the system of coordinates t,ti : x±. 
The lattice used in these computations usually represents a sub-volume of the inter¬ 
action region that expands in the longitudinal direction, as illustrated in the figure 
34. This implies that the lattice spacing in the 2 : coordinate is time dependent. In 
order to be able to resolve the physically relevant scales at the final time Tf of the 
simulation, it is usually necessary to have a larger number of lattice spacings in the 

' One would recover the pathological behavior of the NLO result by linearizing the equation of 
motion for the classical field of initial condition «4i n it + a. 

g Although the discussion here is centered on the energy-momentum tensor, the same resummation 
can be applied to any inclusive quantity. 
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Fig. 34. Lattice setup for numerical implementations of the classical statistical method. 


longitudinal direction. 


10.3. CSA in quantum mechanics 

The Classical Statistical Approximation has an analogue in ordinary quantum me¬ 
chanics, which is helpful to understand what approximation is being made when one 
uses it. The starting point is the transformation that goes from the von Neumann 
equation (36) for the density operator to the Moyal equation (38) for its Wigner 
transform. The Moyal equation is still an exact representation of the quantum 
evolution of the system, although it is expressed in terms of objects that depend 
on classical phase-space variables. However, we have seen that if one expands the 
operator in the right hand side of eq. (38) to lowest order in h , one recovers the 
classical Poisson bracket. In this approximation, the evolution of the system is thus 
purely classical. Note however that the initial condition can remain fully quantum 
in this approximation. For instance, if the system is initially in a pure quantum 
state | 0), its initial Wigner distribution would be 

W 0 {x,p) = J d 3 r e lpr (x + ^ \^)(ip\x - = J d 3 r e lp r ip*{x + ^)ip{x - ^) . 

(63) 

If the initial state is known, then no approximation is needed in the CSA for the 
initial Wigner distribution. In this analogy, the Gaussian initial Wigner distribution 
that we have obtained in the previous subsection would correspond to starting from 
a coherent state, i.e. a state whose Wigner distribution is a Gaussian of width h 
centered on some classical configuration. 


10.4. CSA from the path integral 

It is also instructive to see how the CSA can be derived from the path integral 
formulation of quantum field theory. Since we aim at calculating the expectation 
value of observables in a system initialized in a known state, we must use the 




























September 7, 2015 1:7 


World Scientific Review Volume - 9.75in x 6.5in 


QGP.5 page 45 


45 


Schwinger-Keldysh formalism 

(O) = J [D^{x)\ J [D<j>±(x)} e «W+HW-] 0{(j)) ? (64) 

where the field is doubled into a + component (that represents the evolution in the 
amplitude) and a — component (that describes the conjugate amplitude). In this 
formulation, the initial state of the system is represented by a density operator p 0 , 
whose matrix elements Po(^+\ V'-) control the distribution of the initial values of 
the fields <f>±. The second integral is restricted to fields whose boundary value at 
the initial time is 


(f>±(t 0 ,x) = (p±\x) . (65) 

Inclusive observables do not put any constraint on the final state, and therefore the 
fields have no specific boundary condition at t = +oo in the above path integral, 
except <j> + = 0_ when the measurement is done. 

From eq. (64), one should define new fields as the sum and difference of 0 + and 

0 -, 


^ = 0 + , 0i = 0+ - 0_ . (66) 

The difference 5[0+] — is obviously odd in the field 0i, and it is trivial to 

verify that the term linear in 0 1 comes as a prefactor of the classical equation of 
motion for 02. In addition, in an interacting theory, there are terms that are cubic 
in 0i, 

S'[0_)_] — <S'[0_] = 0i • + terms cubic in 0i . (67) 

002 

We are seeking an approximation that applies in the regime of strong fields, e.g. 
when the fields are excited by a large external source like in the CGC framework, 
which implies that 02 is large. Since 0 + and 0_ are the fields in the amplitude 
and conjugate amplitude respectively, their difference is a quantum effect whose 
magnitude is suppressed by H. Therefore, in such a situation, we have 0i <C 02, 
and it is natural to neglect the cubic term in 0i in the action. After doing this, 
the field 0i becomes a Lagrange multiplier for the classical equation of motion for 
02- The evolution of 0 2 is now deterministic, and the only remaining fluctuations 
are those inherited from the average over the initial density matrix po(<p+\ In 

the Hamiltonian formulation of the classical equation of motion for 02, the average 
over becomes an average over an d its conjugate momentum , with 

a distribution obtained as the Wigner transform of po ■ 

i r v-iio / (*) i x (i) X\ 

J x Pom + 2’^2 - • 


Wirf.nW] = J [D X ] e‘ 


( 68 ) 
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10.5. CSA in perturbation theory 


The path integral derivation of the classical statistical approximation also clarifies 
what this approximation amounts to in perturbation theory. This knowledge will be 
useful later when we discuss the non-renormalizability of this approximation. Let 
us discuss this in the simple framework of a scalar field theory with a </> 4 interaction 
term. 

The first step is to express the transformation from the <f>± fields to the 4>i,i fields 
as a “rotation”, in order to obtain the diagrammatic rules in this new basis. 151-155 
From the propagators in the ± basis, 


G° ++ (P) = 


p 2 — m z + ie 


G°—(p) = 


—i 


p 2 — m? — ie 


G+_(p) = 27 Td(—p°)S(p 2 — m 2 ) , G°_ + (p) = 2tt9(p°)S(p 2 — m 2 ) 


(69) 


one can define a set of new propagators by a linear transformation on the two indices 

g° q/3 = £ n ae n^G^,, (70) 

e,e'=± 

with the following transformation matrix 

= (l /2 I/ 2 ) ■ (?1) 

The free rotated propagators are 

= (Jo G°) ’ (72) 

where we define 

G° r = G° ++ - G° + _ , G° a = G+ + - G°_ + ,G° s = l (G° _ + G°_ + ) . (73) 

(The subscripts R, A and S mean respectively for retarded , advanced and sym¬ 

metric.) Note that the symmetric propagator depends on the initial state of the 
system, while the retarded and advanced ones are independent of this initial data 
(they only reflect how modes propagate in the theory under consideration). Under 
the rotation of eq. (71), the vertices are transformed into 


rim — Tii 22 — 1^2222 — 0 

ri222 = — ig 2 , Tni 2 = — ig 2 / 4 :. (74) 

(The ones not listed explicitly here should be obtained by permutations of the 
indices.) The CSA simply amounts to setting to zero the vertex Tiii 2 (and its 
permutations) wherever it appears in the diagrammatic expansion, while all the 
other Feynman rules are unmodified. 
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11. Applications of the CSA to heavy ion collisions 

The classical statistical approximation has recently been applied to heavy ion colli¬ 
sions in two sets of works, that mostly differ by the nature of the initial conditions 
that were used. 

11.1. CGC Initial conditions 

In a strict application of the CSA to the description of heavy ion collisions in the 
Color Glass Condensate framework, the Gaussian ensemble of initial conditions 
arises from the exponentiation of the 1-loop result. In other words, the CGC de¬ 
scription of heavy ion collisions is an initial value problem, and the state of the 
system at r = 0 + is therefore prescribed uniquely from the fact that the system 
was in the vacuum state at x° = —oo. The initial Gaussian ensemble of fields is 
characterized by the following mean values and variance (written here in a very 
sketchy way, without color and Lorentz indices): 


( A “) = -Ko 



V,V“S^ - V„V' + ig = 0 

lim ak(x) = e lk ' x . 

x°—>— OO 


(75) 


The mean value is already known from the LO calculation , 36 given in eqs. (23). 
From these equations, one sees that the determination of the variance T 2 amounts 


to solving the linearized Yang-Mills equations for all the mode functions (this 
formula comes from the derivation of the NLO contribution ). 157 As long as we stay 


in a regime where these perturbations are not yet enhanced by instabilities, the 
variance is small compared to the mean value squared, and the Gaussian distribution 
of the initial fields is a narrow distribution centered on the LO color fields (see the 
right figure 35). The mode functions describe how plane waves with a well defined 
momentum, color and polarization at x° = —00 are distorted while they propagate 
over the LO color background field . The linearized Yang-Mills equations must 
be solved for each mode function, up to the time tq (in the forward light-cone) at 
which the numerical simulation will start (left figure 35). Explicit formulas for these 
mode functions were derived in ref., 15 ' in the Fock-Schwinger gauge A r = 0 used in 
these computations. For given quantum numbers : v (the Fourier conjugate of the 
rapidity 77 ), fej_,A,c, the gauge potential and the associated electrical fields read 


a i = p +l + p 



e i = -ivfp^-p-A 


e n = -v l (p +i - /r‘) , 


(76) 
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Fig. 35. CGC spectrum of initial fluctuations at Q s tq <C 1. 


where 


/ / 2 \ ^ i 0 

e i P± x ± £j i{p± + k±) / j (gij _ 2 P -^)P X 

P± 

(77) 

and 

= e-^T^)e^Ul(x ± ) j U 2 (p ± + k ± ) (^’-2^)4 . 

P± 

(78) 

U\ 2 are the Wilson lines that describe the color source content of the two colliding 
nuclei (they are defined in terms of p\ 2 in eq. (24)). 

These mode functions have been used in ref. 158 in order to compute the time 
evolution of the components of the energy-momentum tensor from CGC initial 
conditions. The main difficulty and source of uncertainty in these calculation is 
the subtraction of the ultraviolet divergences. Firstly, since the energy-momentum 
tensor is a dimension four operator, it picks up a pure vacuum contribution that 
behaves as the fourth power of the inverse lattice spacing. This contribution can 
be removed by subtracting the result of a second computation done with the same 
lattice parameters, but without the colliding nuclei. A large statistics in the Monte- 
Carlo evaluation of the average over the initial fields is required for this subtraction. 
A second kind of ultraviolet contribution, that depends on the background field, was 
also observed to affect the energy density and the pressure 11 , but not the transverse 
pressure. Since at this level of approximation, the energy momentum tensor is 
traceless: 

_ e = P L +2P T , (79) 

h This contribution is possibly due to the non-renormalizability of the CSA - see the subsection 
12 . 1 . 
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and obeys Bjorken’s law due to energy-momentum conservation, 

d T (re) + P L = 0, (80) 

the only possible form of an ultraviolet sensitive contribution that does not affect 
P T is a term const x r~ 2 that affects equally e and P L . This is indeed what 
was observed. Lacking a more precise understanding of this term, the constant was 
fitted in order to subtract this term from e and P L . After this subtraction, the early 
time behavior of the ratios P L /e and P T /e closely follows the LO up to Q s t ~ 1, 
which is indeed expected since this is before the unstable modes may have affected 
the time evolution. For a moderately large coupling g = 0.5, it was found that 


t [fm/c] 

0.01 0.1 12 3 4 



0.1 1.0 10.0 20.0 30.0 40.0 

Qs i 


Fig. 36. Behavior of P,/e and P. r /e in the CSA with CGC initial conditions. The upper time 
scale, in fm/c , is based on a saturation momentum Q s = 2 GeV. 


the longitudinal pressure increases significantly compared to its value in the LO 
calculation (see the figure 36), and now becomes a sizable fraction of the transverse 
pressure. 


11.2. Particle-like initial conditions 

Alternatively, one may depart somewhat from the CGC framework and use the 
same approximation scheme with a different Gaussian ensemble of initial conditions. 
The simplest model of initial conditions one may consider is an ensemble of fields 
that describe a distribution of free particles. 159-162 In this case, the center of the 
Gaussian is A M = 0, and the variance is constructed from the in-vacuum mode 
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functions (derived in 163 for gluons in the Fock-Schwinger gauge) as follows : 

(A^) = 0 , T 2 (u,v) = J fo(k) a k (u)a* k (v) , a k (x) = e lk ' x , (81) 

modes k 

where fo(k) is the initial momentum distribution of these particles. It is plausible 
that the CGC initial state at r 0 = 0 + may evolve after some time into an incoherent 
distribution of gluons such as the one described by eqs. (81), but this certainly 
deserves a thorough study. 



Fig. 37. Illustration of the difference between the initial conditions (75) and (81) for the Wigner 
distribution of the initial fields. 

The difference and complementarity between the initial conditions in eqs. (75) 
and (81) becomes more transparent if one recalls the symmetric 2-point Green’s 
function in the Schwinger-Keldysh formalism in the presence of a bath of particles, 

G s (k) = 2n^+f 0 (k))5(k 2 ) . (82) 

In this formula, the term in / 0 (fe) represents an initial particle distribution, while 
the term in 1/2 corresponds to pure vacuum quantum fluctuations. It thus be¬ 
comes clear that the CGC initial conditions at r = 0 + can be viewed as vacuum 
fluctuations that are somewhat altered by the presence of the LO color background 
field. It is the minimal modification that quantum mechanics can bring to a clas¬ 
sical state, promoting it to a (pure) coherent state, while the particle-like initial 
conditions defined by eqs. (81) are a mixed state from the point of view of quantum 
mechanics. 

This type of particle-like initial condition was used in refs . 159-162 in order to 
study the evolution of a longitudinally system, in Yang-Mills theory and also in the 
t jA scalar held theory . 162 Note that, if one chooses an initial distribution propor¬ 
tional to the inverse coupling <7 2 , then one can completely scale out the coupling 
from the calculation by an appropriate rescaling of the fields. In these works, the 
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Fig. 38. Behavior of P L /P T in the CSA with particle-like initial conditions in Yang-Mills theory, 
for various anisotropies (£o) and densities (no) of the initial particle distribution. From Ref. 159 

shape of the initial distribution was controlled by two parameters: n o (that con¬ 
trols the overall normalization) and £o (that controls the anisotropy). Some results 
regarding the behavior of the ratio P L /P T are shown in the figure 38. Despite sub¬ 
stantial variations of these parameters, the system was always observed to reach a 
scaling regime in which on has the following approximate behaviors, 

f(p) ~ r~ 2/3 , p ± ~ t° , p z ~ r” 1/3 , ^ ~ t ~ 2/3 . (83) 

Note that these scaling laws are not the free streaming ones (where p z ~ r _1 , 
f(p) ~ t° and P L /P T ~ r -2 ) which means that this system interacts significantly 
but not strongly enough in order to overcome the expansion, at least in the classical 
approximation. The scaling behavior of the longitudinal momentum p z ~ r^ 1 * / 3 
can be understood semi-analytically 112 if the elastic scattering rate is dominated 
by small angle scatterings. Quite surprisingly, the same scaling behavior of p z was 
also observed 162 in this approximation in the scalar </> 4 theory, despite the fact that 
the leading order scattering cross-section in this theory is dominated by large angle 
scatterings 1 . 

In this approximation, this scaling behavior of eq. (83) will continue forever, and 
the pressure tensor becomes more and more anisotropic over time. However, since 
this is a classical approximation, it is known to break at least when the occupation 
number f(p) becomes of order one. This is expected to happen when Q s t ~ 
g~ 3 , which corresponds to the end of the first stage of the bottom-up scenario of 
ref. 112 If this scenario illustrated in the figure 39- is confirmed, the isotropization 

1 See the section 12.3 for a discussion of the possible interplay between the classical approximation 

and large angle scatterings in anisotropic systems. 
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Fig. 39. Thermalization/isotropization scenario that emerges from CSA studies with particle-like 
initial conditions. From Ref. 159 

of the system would only begin after this time, while the occupation number in 
the system is no longer large. The quantum corrections would eventually bring 
the energy-momentum tensor to a quasi-isotropic shape compatible with nearly 
ideal hydrodynamics. However, what would be more difficult to understand in this 
scenario is why the shear viscosity to entropy ratio is small. Indeed, in a weak 
coupling system with occupation numbers of order 1 or below, this ratio is expected 
to be parametrically large g/s ~ g~ A . 

12. Limitations of the classical statistical approximation 

The main appeal of the CSA is the fact that it is quite straightforward to implement, 
(unlike other schemes like the two-particle irreducible approach, that has never been 
applied to Yang-Mills theory so far), while at the same time staying very close to the 
dynamics of the gauge fields so that it is a natural extension of the CGC framework. 
However, it is not without problems and limitations, that we discuss in this section. 

12.1. Ultraviolet divergences and non-renormalizability 

A very important difference between the initial conditions (75) and (81) is the 
momentum dependence of the spectrum of field fluctuations. In eqs. (81) the large 
momentum behavior of the spectrum is controlled by the initial particle distribution 
/o(fc), and therefore it does not extend to infinity if one chooses a /o(fe) that has a 
compact support. In contrast, the vacuum fluctuations that are the source of the 
field fluctuations in the CGC initial conditions have a spectrum which is flat up to 
k = oo. 

This difference in the spectra at large momentum leads to very different de- 
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pendences on the ultraviolet cutoff (i.e. the inverse lattice spacing) in numerical 
implementations of the CSA using these two types of initial conditions. This issue 
can be phrased as follows: starting from a renormalizable quantum field theory (e.g. 
Yang-Mills theory, or a (j) A scalar field theory), does the CSA preserve its renormal- 
izability? It turns out that the answer to this question depends on the spectrum 
of the initial field fluctuations. This can be studied perturbatively 155 by using the 
retarded-advanced basis, and by using the fact that the CSA amounts to dropping 
the vertex that has 3 indices of type 1. 

From studies performed in the context of quantum field theory at finite tem¬ 
perature, 154,164-16 ' it has been known for a long time that a particle-like spectrum 
that falls at least as fast as fc -1 leads to a super-renormalizable approximation. 
In this case, it is sufficient to perform a finite number of subtractions in order to 
make predictions that do not depend on the ultraviolet cutoff. This applies directly 
to initial conditions of type (81), provided that the initial distribution /o(fc) falls 
quickly enough with momentum. 

The situation is quite different with the vacuum-like CGC initial conditions, 
because they have a flat spectrum of fluctuations. By using the perturbative CSA 
described in the subsection 10.5, it has been shown in Ref. 150 that the CSA is 
a non renormalizable approximation of the underlying quantum field theory when 
this type of initial condition is used. For instance, for the self-energy E 12 at two 
loops (in a </> 4 scalar field theory), the CSA gives an ultraviolet divergent result 


Im 



1024tt 3 



(84) 


Since this divergence occurs in the imaginary part of a correlator, it cannot be re¬ 
moved by a counterterm added to the action (otherwise that would make the action 
non Hermitian). The consequence of this is that one cannot take the continuum 
limit in computations based on the CSA with vacuum-like initial conditions. This 
2-loop divergence in a self-energy has a counterpart in the classical approximation 
of the Boltzmann equation, as we shall see in the next subsection. 

Non-renormalizable graphs can even be found in some 1-loop 4-point functions, 
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such as the following Fn 22 function : 


-iT 


1 loop 
1122 



S channel 




(85) 


Here, we have shown all the vertex assignments that would appear in the retarded- 
advanced basis. Calculating and adding up all these contributions would lead to 
an ultraviolet finite result, as is expected in a renormalizable theory since the bare 
Lagrangian has no 4> 2 </>| operator. Note however that some of these graphs contain 
a 1112 vertex, and would therefore be discarded in the CSA. In this approximation, 
we only have 



the result of which is given by 

„4 ■ 

sign (f) 


• r-p 1 1 loo P 
— * |1 H22j CSA 


9 

64-7T 


■ sign(u) + 2 A„ 


0(-t) 

\Pi +P 3 | 


with 


( 86 ) 


g(-«) y 

\pi+Pi\)\ 

(87) 


t = (p 1 +p 3 ) 2 , n=(pi+p 4 ) 2 , (88) 

the standard Mandelstam variables and where A uv is the ultraviolet cutoff on 3- 
momentum. Despite the zero superficial degree of divergence of these graphs, they 
contain a linear divergence. Moreover, the coefficient of the divergent terms is non¬ 
polynomial in the momenta, implying that it is a non-local ultraviolet divergence. 
Such non-renormalizable contributions will appear at 2-loops and beyond in the 
expectation value of inclusive quantities like the energy-momentum tensor. Because 
of them, CSA calculations performed with initial conditions that contain vacuum 
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fluctuations should be performed with an ultraviolet cutoff (i.e. the inverse lattice 
spacing) which is not too large compared to the physical scales. 


12.2. CSA in kinetic theory 

Assessing the cutoff dependence within the CSA itself is very costly, because it 
requires repeating the same calculation several times with smaller and smaller lat¬ 
tice spacings. 161 A much cheaper way of understanding the interplay between the 
classical approximation and the dependence on the ultraviolet cutoff is to consider 
the same approximation at the level of kinetic theory 1 . Let us start with the Boltz¬ 
mann equation, with the collision term expressed in terms of self-energies in the 
Schwinger-Keldysli formalism : 


[dt + v P ■ V] f( P ) = 2 ^- [/(P)S_+(P) - (1 + /(p))£ + _(P)] . (89) 

CAf] 


In order to perform in kinetic theory the same approximation as in the CSA, one 
should first rewrite the collision term in the retarded-advanced basis : 


Cplf] 




Sn(P)+(/(p )+2 ) (S 2 i(P) - Si 2 (P)) 


(90) 


(£n is imaginary, as well as £ 21 (P) - £12(P ))■ At this point, we just need to 
calculate the self-energies that appear in the right hand side (at 2-loops if we want 
2 —>■ 2 collisions) by neglecting the 1112 vertex. The two versions of the CSA, 
with or without vacuum fluctuations, simply correspond to keeping the 1/2 in the 
factor 1/2 + f(p) that appears in the propagator G 22 or not. 174-176 It is easy to 
check that the classical approximation with no vacuum fluctuations amounts to 
keeping only the terms that are cubic in the particle distribution, while the classical 
approximation with vacuum fluctuations leads to the following collision term 


T~ (2n) 4 S(P + K — P' — K') 

4cJp Jk Jp' Jk' 

x[(f(p') + l)(f(k') + ±)(l + f(p)+f(k)) 


(91) 


that has the same cubic and quadratic terms as the exact collision term, but also 
some spurious terms that are linear in the particle distribution. 

J Here, we employ kinetic theory as a way to assess some formal aspects of the underlying quantum 
field theory, like its dependence on the ultraviolet cutoff. A number of works 108-111,168-173 have 
also used kinetic theory as a tool for studying thermalization in models of heavy ion collisions. 
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The fixed points of these “classical” Boltzmann equations are 

T 

without vacuum fluctuations : f(p ) = - 

^p 

T 1 

with vacuum fluctuations : f(p) = -- , (92) 

c Op jU 2 

(LOp = Vp* + m 2 ) i.e. respectively the first and the first two terms in the expansion 
of a Bose-Einstein distribution around low energy. 

The parameters T and p that appear in the asymptotic distributions can be 
determined from conservation laws. In the Boltzmann equation with only elastic 
scatterings, both the energy and the particle number are conserved. For the CSA 
with vacuum fluctuations, these conservation laws lead to 


n = n c + 


e = n c m + 


d 3 p 

( 2^)3 

f d 3 p 


T 


LO p - fi 


(27r) 


3 P 


T 


LOp - p 


(93) 


In this exercise, we allow the formation of a Bose-Einstein condensate with a particle 
density n c , in case one wishes to consider highly populated initial conditions. Note 
that in any case, there are only two unknowns in these equations: T and p if 
there is no condensation, or T and n c if there is condensation (in which case p = 
m). From these two equations, we can determine the two unknown parameters. 
However, since the integrals are ultraviolet divergent, it is necessary to introduce 
a cutoff A on |p|, which turns T, p and n c into A-dependent quantities. 156 The 


m = 0.5 0 £ = Q 4 n = 0.75 £ / m [ Classical + 1/2 ] 



Fig. 40. Evolution of T, fi and n c as a function of the ultraviolet cutoff. The points reproduce 
the values listed in the figure 10 of Ref., 161 obtained with a classical statistical field simulation. 

solution of eqs. (93) is shown in the figure 40 (Q is a physical momentum scale 
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that characterizes the initial momentum distribution of the particles), where we 
have also superimposed results from ref . 161 obtained by lattice classical statistical 
simulations. These curves show a very strong dependence on the ultraviolet cutoff 
when it becomes much larger than the physical scales, as expected given the fact 
that the CSA with vacuum fluctuations is not renormalizable. On the other hand, 
there is a region where the cutoff is a few times the physical scale, and where the 
parameters that characterize the asymptotic distribution are rather insensitive to 
the cutoff. Therefore, simulations using the CSA with vacuum fluctuations as initial 
conditions should preferably be performed with an ultraviolet cutoff chosen in that 
range, in order to minimize the sensitivity of the results on the cutoff. 

12.3. Quantum, corrections in anisotropic systems 

Note that this non-renormalizability problem arises only when one uses the CSA 
in conjunction with a spectrum of initial conditions that represents vacuum fluc¬ 
tuations, as in eqs. (75). When using particle-like initial conditions, such as those 
described in eqs. (81), the CSA is ultraviolet finite provided that the initial particle 
distribution falls faster than 1/k. Because of this, if one forgets that this type of 
initial condition is not derived directly from the CGC, this implementation of the 
CSA may seem better since one does not need to worry about the dependence on 
the ultraviolet cutoff. 

However, this choice of initial conditions in the CSA leads to a different kind 
of problem when used in studies of isotropization in heavy ion collisions, for the 
following reason. For the purpose of this argument, let us again reason in terms of 
the Boltzmann equation. For 2 — > 2 collisions, it reads 

dth ~ g 4 [ ■■■ [/i/ 2 (/ 3 + h) - hhih + h)] 

J 124 

+g 4 f ■■■ [/ 1/2 - hh] ■ (94) 

J 124 

(Here we are tracking the distribution of particles of momentum p 3 - see the figure 
41.) In this equation, we have separated the purely classical terms (in / 3 ) from the 
subleading / 2 terms. When only particle-like fluctuations are included, the CSA 
applied to the Boltzmann equation keeps only the / 3 terms, and neglects all the 
other terms. In contrast, the CSA where vacuum fluctuations are included has the 
/ 3 and / 2 terms, and also some unphysical terms that are linear in / (see Eq. (91)). 

Consider now a situation where the particle distribution is strongly anisotropic, 
with a support in p. which is squeezed compared to the support in p±. Starting 
from a generic initial condition of this form, the unapproximated Boltzmann equa¬ 
tion will in general lead to isotropization because two purely transverse particles 
can be scattered outside of the transverse plane k . But does this still happen in ap¬ 
proximations of the Boltzmann equation? In the figure 41, one has = —p% 7 ^ 0, 


k This is obvious in the ct' 1 scalar field theory, where the leading term in the two body cross-section 
is point-like and where scatterings at large angle are dominant. This is also the case in the CGC, 
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Fig. 41. 2 —>• 2 scattering contributing to isotropization. 


and given our assumption about the support of the particle distribution, this means 
that fs = fi S3 0. In the collision term of the Boltzmann equation, a number of 
terms therefore vanish : 


dth ~g 4 J ■■■ [fiM(3+jA) -hhifi + h)\ 
124 0 0 

+g 4 [ • • • [/ 1/2 - 73 / 4 ] ■ 

J 124 '-v'-'' 


(95) 


In particular, all the cubic terms are zero. The problem is that these are the only 
terms that are kept in the CSA with no vacuum fluctuations. The only non-zero 
term is the term in / 1 / 2 , which would be present in the CSA only if vacuum initial 
fluctuations are present. From this discussion, the particle-like initial conditions in 
the CSA, despite their appeal since they lead to UV finite results, may be inappropri¬ 
ate because they lead to missing the most important contribution to isotropization. 
In other words, in an anisotropic system, the classical approximation could break 
down long before the occupation number becomes of order 1. 

This may also have an incidence on the behavior at intermediate times. With 
only the / 3 terms, the ratio P L /P T decreases like t _2//3 at late times, and therefore 
the system never isotropizes in this approximation. On the other hand, the system 
is expected to isotropize eventually with the complete Boltzmann equation (/ 3 and 
/ 2 terms), with a ratio P L /P T ~ r° at late times. If the f 2 terms are truly negligible 
over some extended period of time, then the full Boltzmann equation should lead 
to the red curve in the figure 42, in which the system spends some time stuck into a 


thanks to screening. Indeed, since the occupation number is of order g 2 , the Debye mass is of 
order Q s . Since this is comparable to the typical momentum of the gluons, large angle scatterings 
should be important at early times. 
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, dlog(P L /P T ) 

dT 



Fig. 42. Behavior of the logarithmic derivative of the ratio P L /P T . Dotted curve: classical 
approximation where one keeps only the / 3 terms. Red curve: full collision term, if there exists a 
“classical attractor”. Blue curve: full collision term, if the / 2 terms prevent the classical attractor. 

“classical attractor” before eventually leaving it in order to isotropize. In contrast, 
if the fact that the f 2 terms are always dominant in the tail plays an important 
role, then one may instead get the blue curve. The Boltzmann equation offers an 
interesting playground in order to test these possibilities, since it can be solved with 
and without the / 2 terms at comparable computational costs. 

13. Summary 

A lot of progress has been made in the past 10 years in QCD-based studies of 
the early stages of heavy ion collisions. The most promising framework for these 
studies is the Color Glass Condensate, which allows one to systematically include 
the non-linear effects that prevail when the gluon occupation numbers are large. 

The description of the relevant degrees of freedom in the wavefunction of the 
incoming nuclei and how they evolve as one varies the energy of the collision is 
entering a very mature stage, since one now knows this evolution at next-to-leading 
log accuracy. A lot more work will be needed before these NLL evolution equations 
can be implemented in phenomenological studies of heavy ion collisions, but this is 
a very important step towards more quantitative results. 

Regarding the collision itself, i.e. how two objects described by means of the 
CGC interact while they collide, one has now a much clearer understanding of how 
observables can be expressed at leading order in terms of a classical solution of the 
Yang-Mills equations (and how the inclusive nature of an observable translates into 
retarded boundary conditions for this classical solution), and how to calculate their 
next-to-leading order corrections in terms of linearized perturbations around this 
classical solution. 

The recent years have also witnessed a renewed interest in the question of the 
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isotropization and thermalization of the system produced in heavy ion collisions. 
This issue has indeed been made more pressing by the many phenomenological 
successes of hydrodynamical models in describing the expansion of the fireball as a 
nearly perfect fluid, something which is at odds with the LO CGC results. From 
recent works, there appears to be two tools of choice for these studies: kinetic theory 
(in the small scattering angle approximation, or the more sophisticated effective 
kinetic theory of ref. 105 ) and the classical statistical approximation, which is more 
directly related to the CGC framework. Although this problem seems to have now 
reached a satisfactory level of understanding -both qualitatively and quantitatively- 
in the simpler case of a system confined in a fixed volume, the situation is at 
the moment somewhat unclear in the more realistic situation of a longitudinally 
expanding system, for which two different results have been obtained in the classical 
statistical approximation with different initial conditions. These initial conditions 
mostly differ in whether some vacuum quantum fluctuations are included or not, 
and a prime question to clarify in the near future will be to understand to what 
extent quantum fluctuations play a role in isotropization and thermalization. 
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